Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Le problème d’apprentissage

Le chapitre précédent a présenté les méthodes non paramétriques, comme les k plus proches voisins, qui conservent les données d’entraînement et les consultent au moment de la prédiction. Cette approche est intuitive, mais elle a un coût: les données doivent rester en mémoire, et chaque prédiction requiert de parcourir l’ensemble d’entraînement. Ce chapitre développe une approche différente: plutôt que de garder les données, nous cherchons à les résumer dans un ensemble fixe de paramètres. L’apprentissage devient alors un problème d’optimisation.

Apprentissage supervisé

Une ingénieure automobile mesure la distance de freinage d’un véhicule à différentes vitesses. Ses données ressemblent à ceci:

Vitesse (mph)Distance (ft)
42
74
1220
1856
2493

Elle veut prédire la distance de freinage à 30 mph sans faire l’essai. Pour cela, elle cherche une fonction ff telle que f(vitesse)distancef(\text{vitesse}) \approx \text{distance} sur ses observations. Si la fonction capture la relation sous-jacente, elle devrait donner une prédiction raisonnable pour des vitesses non mesurées.

Source
import numpy as np
import matplotlib.pyplot as plt

# Configuration pour des figures haute résolution
%config InlineBackend.figure_format = 'retina'

# Données de freinage (Ezekiel, 1930): vitesse (mph) vs distance d'arrêt (ft)
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

plt.figure(figsize=(6, 4))
plt.scatter(speed, dist, alpha=0.7, label='Observations')

# Fit quadratic
coeffs = np.polyfit(speed, dist, 2)
speed_grid = np.linspace(0, 30, 100)
dist_pred = np.polyval(coeffs, speed_grid)
plt.plot(speed_grid, dist_pred, 'k--', alpha=0.6, label='Fonction ajustée')

# Prediction at 30 mph
pred_30 = np.polyval(coeffs, 30)
plt.scatter([30], [pred_30], marker='x', s=80, color='C1', zorder=5, label=f'Prédiction à 30 mph: {pred_30:.0f} ft')

plt.xlabel('Vitesse (mph)')
plt.ylabel('Distance de freinage (ft)')
plt.legend()
plt.tight_layout()
<Figure size 600x400 with 1 Axes>

Ce processus est l’ajustement de courbe (curve fitting). Nous avons des paires (entrée, sortie), nous ajustons une fonction, et nous utilisons cette fonction pour prédire. L’apprentissage supervisé généralise cette idée: les entrées peuvent être des vecteurs de dimension quelconque, les sorties peuvent être continues ou discrètes, et les fonctions candidates peuvent être bien plus complexes qu’un polynôme.

Formellement, nous disposons d’un jeu de données D={(xi,yi)}i=1N\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^N composé de NN paires, où chaque xiX\mathbf{x}_i \in \mathcal{X} est une entrée et yiYy_i \in \mathcal{Y} est la sortie correspondante. L’objectif est de trouver une fonction f:XYf: \mathcal{X} \to \mathcal{Y} qui approxime bien la relation entre entrées et sorties, y compris pour des exemples que nous n’avons pas encore observés.

Dans de nombreuses applications, les entrées sont des vecteurs de caractéristiques. Chaque exemple xiRd\mathbf{x}_i \in \mathbb{R}^d est un vecteur de dimension dd, où chaque composante représente une mesure ou un attribut. Pour prédire le prix d’une maison, les entrées pourraient être la superficie, le nombre de chambres et l’âge du bâtiment. Pour filtrer les pourriels, les entrées pourraient être des fréquences de mots. Pour diagnostiquer une maladie, les entrées pourraient être des résultats d’analyses sanguines.

Lorsque la sortie est une valeur continue, nous parlons de régression: f:RdRf: \mathbb{R}^d \to \mathbb{R} pour une sortie scalaire, ou f:RdRpf: \mathbb{R}^d \to \mathbb{R}^p pour une sortie vectorielle. La distance de freinage, le prix d’une maison, la concentration d’un médicament dans le sang sont des exemples de régression.

Lorsque la sortie est une catégorie parmi un ensemble fini, nous parlons de classification. Pour la classification binaire, f:Rd{0,1}f: \mathbb{R}^d \to \{0, 1\}. Pour la classification multiclasse avec mm catégories, f:Rd{0,,m1}f: \mathbb{R}^d \to \{0, \ldots, m-1\}. Déterminer si un courriel est un pourriel, diagnostiquer une maladie, ou reconnaître un chiffre manuscrit sont des exemples de classification.

Mesurer l’erreur

Pour choisir entre deux fonctions candidates, nous avons besoin d’un critère qui quantifie la qualité des prédictions. Une fonction de perte :Y×YR+\ell: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}_+ mesure l’écart entre une prédiction y^\hat{y} et la vraie valeur yy. Une perte de zéro indique une prédiction parfaite; plus la perte est grande, plus l’erreur est importante.

Pour la régression, nous utilisons généralement la perte quadratique:

2(y,y^)=(yy^)2\ell_2(y, \hat{y}) = (y - \hat{y})^2

Cette perte pénalise les grandes erreurs de manière quadratique. Une erreur de 2 coûte quatre fois plus qu’une erreur de 1.

Reprenons les données de freinage. Supposons que notre fonction prédise 50 ft pour une vitesse où la vraie distance est 56 ft. La perte quadratique est (5650)2=36(56 - 50)^2 = 36. Si elle prédit 70 ft, la perte est (5670)2=196(56 - 70)^2 = 196. La perte quadratique pénalise sévèrement les grandes erreurs.

Source
import numpy as np
import matplotlib.pyplot as plt

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

coeffs = np.polyfit(speed, dist, 2)
predictions = np.polyval(coeffs, speed)
residuals = dist - predictions

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Left: predictions vs observations
ax = axes[0]
ax.scatter(speed, dist, alpha=0.7, label='Observations')
speed_grid = np.linspace(4, 25, 100)
ax.plot(speed_grid, np.polyval(coeffs, speed_grid), 'k--', alpha=0.6, label='Prédictions')
for i in range(len(speed)):
    ax.plot([speed[i], speed[i]], [dist[i], predictions[i]], 'C1-', alpha=0.5)
ax.set_xlabel('Vitesse (mph)')
ax.set_ylabel('Distance (ft)')
ax.legend()
ax.set_title('Résidus: écarts entre observations et prédictions')

# Right: histogram of squared residuals
ax = axes[1]
ax.hist(residuals**2, bins=15, edgecolor='black', alpha=0.7)
ax.set_xlabel(r'Perte quadratique $(y - \hat{y})^2$')
ax.set_ylabel('Fréquence')
ax.set_title(f'EQM = {np.mean(residuals**2):.1f}')

plt.tight_layout()
<Figure size 1000x400 with 2 Axes>

Pour la classification, un choix naturel est la perte 0-1:

01(y,y^)=1yy^={0si y=y^1si yy^\ell_{0-1}(y, \hat{y}) = \mathbb{1}_{y \neq \hat{y}} = \begin{cases} 0 & \text{si } y = \hat{y} \\ 1 & \text{si } y \neq \hat{y} \end{cases}

Cette perte compte simplement les erreurs: elle vaut 1 pour une mauvaise prédiction, 0 sinon.

Le choix de la fonction de perte dépend du problème. En diagnostic médical, manquer une maladie grave (faux négatif) peut avoir des conséquences bien plus importantes que de prescrire un test supplémentaire à un patient sain (faux positif). Une perte asymétrique refléterait cette différence. En régression, si les grandes erreurs sont particulièrement problématiques, la perte quadratique est appropriée; si nous voulons être robustes aux valeurs aberrantes, la perte absolue yy^|y - \hat{y}| est préférable.

Le risque

La perte évalue une seule prédiction. Pour évaluer un modèle dans son ensemble, nous voulons mesurer sa performance moyenne sur toutes les données possibles, pas seulement sur les exemples que nous avons observés.

Pourquoi des variables aléatoires?

Une question naturelle se pose: si nous ajustons une fonction déterministe ff à des données, pourquoi avons-nous besoin de variables aléatoires et d’espérances? La fonction obtenue n’est-elle pas simplement une courbe fixe?

La réponse tient en un mot: généralisation. Nous ne nous intéressons pas vraiment à la performance sur les données d’entraînement car ces points sont déjà connus. Ce qui compte, c’est la performance sur des données futures que nous n’avons pas encore observées.

Considérons l’exemple de la distance de freinage. Les 50 mesures dans notre tableau sont un échantillon de toutes les mesures possibles. Si nous retournions sur le terrain et mesurions à nouveau, nous obtiendrions des valeurs légèrement différentes. En effet, le même véhicule à 20 mph ne s’arrête pas exactement à la même distance à chaque essai. Il y a de la variabilité intrinsèque: état de la route, température des freins, réflexes du conducteur.

Cette variabilité est capturée par une distribution de probabilité p(x,y)p(x, y). Nos 50 points sont des tirages de cette distribution. La question fondamentale devient alors:

Notre modèle ff prédira-t-il bien sur de nouveaux tirages de cette même distribution?

La fonction ff elle-même est déterministe une fois entraînée. Mais son évaluation, savoir si elle prédit bien ou mal, dépend de quelles données futures elle rencontrera. Et ces données futures sont incertaines: elles seront tirées de p(x,y)p(x, y), mais nous ne savons pas lesquelles.

Le risque formalise cette idée: c’est la perte moyenne que subira notre modèle ff lorsqu’il sera confronté à des données tirées de p(x,y)p(x, y). C’est une mesure de performance prospective, tournée vers le futur.

Définition formelle

Le risque d’une fonction ff est l’espérance de la perte sur la distribution des données:

R(f)=E(X,Y)p[(Y,f(X))]=(y,f(x))p(x,y)dxdy\mathcal{R}(f) = \mathbb{E}_{(\mathbf{X},Y) \sim p}\left[\ell(Y, f(\mathbf{X}))\right] = \int \ell(y, f(\mathbf{x})) \, p(\mathbf{x}, y) \, d\mathbf{x} \, dy

Décomposons cette formule étape par étape:

  1. E(X,Y)p\mathbb{E}_{(\mathbf{X},Y) \sim p}: L’espérance mathématique signifie “moyenne sur tous les exemples possibles”. La notation (X,Y)p(\mathbf{X},Y) \sim p indique que nous tirons les paires (x,y)(\mathbf{x}, y) selon la distribution p(x,y)p(\mathbf{x}, y) de la nature.

  2. (Y,f(X))\ell(Y, f(\mathbf{X})): Pour chaque exemple aléatoire (X,Y)(\mathbf{X}, Y), nous calculons la perte entre la vraie valeur YY et la prédiction f(X)f(\mathbf{X}) du modèle.

  3. L’intégrale (y,f(x))p(x,y)dxdy\int \ell(y, f(\mathbf{x})) \, p(\mathbf{x}, y) \, d\mathbf{x} \, dy: Cette intégrale calcule une moyenne pondérée. Pour chaque paire possible (x,y)(\mathbf{x}, y), nous multiplions la perte (y,f(x))\ell(y, f(\mathbf{x})) par la probabilité p(x,y)p(\mathbf{x}, y) que cette paire apparaisse dans la nature, puis nous sommons (intégrons) sur toutes les paires possibles.

Exemple concret

Considérons un problème de classification binaire en 2D. Supposons que x[0,1]2\mathbf{x} \in [0, 1]^2 et y{0,1}y \in \{0, 1\}. Pour calculer le risque, nous devrions:

  1. Diviser l’espace [0,1]2[0,1]^2 en une grille fine (par exemple, 1000×10001000 \times 1000 points)

  2. Pour chaque point x\mathbf{x} de la grille, considérer les deux valeurs possibles de yy (0 et 1)

  3. Pour chaque combinaison (x,y)(\mathbf{x}, y), calculer:

    • La probabilité p(x,y)p(\mathbf{x}, y) que cette combinaison apparaisse

    • La perte (y,f(x))\ell(y, f(\mathbf{x})) si notre modèle prédit f(x)f(\mathbf{x})

  4. Faire la somme pondérée: xy{0,1}(y,f(x))p(x,y)\sum_{\mathbf{x}} \sum_{y \in \{0,1\}} \ell(y, f(\mathbf{x})) \cdot p(\mathbf{x}, y)

En pratique, pour un espace continu, cette somme devient une intégrale sur un domaine continu, ce qui est encore plus complexe à calculer.

Visualisons ceci concrètement. La figure suivante montre un problème de classification binaire où chaque classe suit une distribution gaussienne en 2D. Les contours représentent la densité p(xy)p(x|y) pour chaque classe. La ligne pointillée est la frontière de décision d’un classificateur linéaire. Les régions ombrées indiquent où le classificateur fait des erreurs: la région rouge correspond aux points de classe 0 classés comme classe 1, et la région bleue correspond aux points de classe 1 classés comme classe 0.

Source
import numpy as np
import matplotlib.pyplot as plt

# Paramètres du mélange gaussien (classification binaire 2D)
mu0 = np.array([0.0, 0.0])
mu1 = np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])

def gaussian_pdf(x, mu, cov):
    """PDF d'une gaussienne multivariée."""
    d = len(mu)
    diff = x - mu
    cov_inv = np.linalg.inv(cov)
    mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
    return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))

# Create grid for visualization
x_range = np.linspace(-3, 5, 200)
y_range = np.linspace(-3, 4, 200)
X_grid, Y_grid = np.meshgrid(x_range, y_range)
pos = np.dstack([X_grid, Y_grid])

# Compute class-conditional densities
p_x_given_0 = gaussian_pdf(pos, mu0, cov)
p_x_given_1 = gaussian_pdf(pos, mu1, cov)

# Joint densities (with equal priors)
prior = 0.5
p_x_y0 = p_x_given_0 * (1 - prior)
p_x_y1 = p_x_given_1 * prior

# Linear decision boundary (Bayes optimal for equal covariances)
# w^T x + b = 0 where w = Sigma^{-1}(mu1 - mu0)
cov_inv = np.linalg.inv(cov)
w = cov_inv @ (mu1 - mu0)
b = -0.5 * (mu1 @ cov_inv @ mu1 - mu0 @ cov_inv @ mu0)

# Decision boundary: w[0]*x + w[1]*y + b = 0  =>  y = -(w[0]*x + b)/w[1]
x_boundary = np.linspace(-3, 5, 100)
y_boundary = -(w[0] * x_boundary + b) / w[1]

# Classifier prediction: classify as 1 if w^T x + b > 0
predictions = (w[0] * X_grid + w[1] * Y_grid + b) > 0

# Misclassification regions
# Class 0 misclassified as 1: true class is 0, but prediction is 1
misclass_0 = predictions  # region where we predict 1
# Class 1 misclassified as 0: true class is 1, but prediction is 0
misclass_1 = ~predictions  # region where we predict 0

fig, ax = plt.subplots(figsize=(8, 6))

# Plot class-conditional densities as contours
levels = [0.01, 0.05, 0.1, 0.15]
ax.contour(X_grid, Y_grid, p_x_given_0, levels=levels, colors='C0', alpha=0.7)
ax.contour(X_grid, Y_grid, p_x_given_1, levels=levels, colors='C1', alpha=0.7)

# Shade misclassification regions weighted by probability
# Red: class 0 points incorrectly classified as 1
error_region_0 = np.where(misclass_0, p_x_y0, 0)
# Blue: class 1 points incorrectly classified as 0  
error_region_1 = np.where(misclass_1, p_x_y1, 0)

ax.contourf(X_grid, Y_grid, error_region_0, levels=[0.001, 0.01, 0.05, 0.1], 
            colors=['#ff000010', '#ff000030', '#ff000050'], extend='max')
ax.contourf(X_grid, Y_grid, error_region_1, levels=[0.001, 0.01, 0.05, 0.1],
            colors=['#0000ff10', '#0000ff30', '#0000ff50'], extend='max')

# Decision boundary
ax.plot(x_boundary, y_boundary, 'k--', linewidth=2, label='Frontière de décision')

# Class centers
ax.scatter(*mu0, s=100, c='C0', marker='x', linewidths=3, zorder=5, label='Centre classe 0')
ax.scatter(*mu1, s=100, c='C1', marker='x', linewidths=3, zorder=5, label='Centre classe 1')

ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.legend(loc='upper left')
ax.set_aspect('equal')

plt.tight_layout()
<Figure size 800x600 with 1 Axes>

Le risque est l’intégrale de la perte sur tout l’espace, pondérée par p(x,y)p(\mathbf{x}, y). Les régions ombrées contribuent au risque: chaque point dans ces régions est mal classé, et sa contribution dépend de la densité de probabilité à cet endroit. Les régions denses proches de la frontière contribuent le plus au risque.

Pourquoi le risque est important

Le risque mesure ce que nous obtiendrons en moyenne si nous appliquons ff à de nouvelles données tirées de la même distribution. Un modèle avec un faible risque fait de bonnes prédictions en général, pas seulement sur les exemples d’entraînement. C’est exactement ce que nous voulons optimiser: un modèle qui performe bien sur des données jamais vues, pas seulement sur celles qu’il a déjà observées.

Cette quantité est ce que nous voulons minimiser. Le problème fondamental est que nous ne connaissons pas la distribution p(x,y)p(\mathbf{x}, y) de la nature. Nous n’y avons accès qu’indirectement, via un échantillon fini D\mathcal{D}.

Le risque empirique

Puisque le risque est inaccessible, nous l’approximons par une moyenne sur les données disponibles. Le risque empirique est:

R^(f,D)=1Ni=1N(yi,f(xi))\hat{\mathcal{R}}(f, \mathcal{D}) = \frac{1}{N} \sum_{i=1}^{N} \ell(y_i, f(\mathbf{x}_i))

Cette quantité est calculable: c’est la moyenne des pertes sur l’échantillon d’entraînement. Pour la perte 0-1, le risque empirique est le taux d’erreur sur les données d’entraînement. Pour la perte quadratique, c’est l’erreur quadratique moyenne.

Pourquoi le risque est-il inaccessible?

La nécessité d’utiliser le risque empirique découle de deux obstacles fondamentaux, l’un conceptuel et l’autre computationnel.

Obstacle 1: La distribution p(x,y)p(\mathbf{x}, y) est inconnue

La nature possède une distribution p(x,y)p(\mathbf{x}, y) qui génère les données, mais nous ne la connaissons pas. Nous n’observons qu’un échantillon fini D={(xi,yi)}i=1N\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^N tiré de cette distribution.

L’ensemble D\mathcal{D} est une variable aléatoire: si nous répétions l’expérience de collecte de données, nous obtiendrions un échantillon différent. Cette perspective, adoptée notamment dans Murphy (2022), rappelle que nos conclusions dépendent de l’échantillon particulier que nous avons observé. C’est comme si nous regardions quelques gouttes d’eau d’un océan: nous pouvons analyser ces gouttes, mais un autre prélèvement donnerait des gouttes différentes.

Même si nous tentions d’estimer p(x,y)p(\mathbf{x}, y) à partir des données (par exemple, via des techniques d’estimation de densité comme les mélanges de gaussiennes ou les estimateurs à noyau), nous n’obtiendrions qu’une approximation p^(x,y)\hat{p}(\mathbf{x}, y) de la vraie distribution. Cette approximation serait elle-même imparfaite et dépendrait de nos hypothèses sur la forme de la distribution.

La figure suivante illustre ce problème. À gauche, la vraie distribution p(x,y)p(\mathbf{x}, y) que la nature utilise pour générer les données (que nous ne connaissons pas). À droite, un échantillon de N=50N = 50 points tirés de cette distribution (ce que nous observons).

Source
import numpy as np
import matplotlib.pyplot as plt

# Paramètres du mélange gaussien
mu0, mu1 = np.array([0.0, 0.0]), np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])

def gaussian_pdf(x, mu, cov):
    d = len(mu)
    diff = x - mu
    cov_inv = np.linalg.inv(cov)
    mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
    return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))

# Générer un échantillon
rng = np.random.default_rng(42)
n = 50
X0 = rng.multivariate_normal(mu0, cov, n // 2)
X1 = rng.multivariate_normal(mu1, cov, n // 2)
X = np.vstack([X0, X1])
y = np.concatenate([np.zeros(n // 2), np.ones(n // 2)])

# Grille pour visualisation
x_range = np.linspace(-3, 5, 150)
y_range = np.linspace(-3, 4, 150)
X_grid, Y_grid = np.meshgrid(x_range, y_range)
pos = np.dstack([X_grid, Y_grid])

p_x_given_0 = gaussian_pdf(pos, mu0, cov)
p_x_given_1 = gaussian_pdf(pos, mu1, cov)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Left: True distribution (what nature knows)
ax = axes[0]
ax.contourf(X_grid, Y_grid, p_x_given_0, levels=15, cmap='Blues', alpha=0.6)
ax.contourf(X_grid, Y_grid, p_x_given_1, levels=15, cmap='Oranges', alpha=0.6)
ax.contour(X_grid, Y_grid, p_x_given_0, levels=5, colors='C0', alpha=0.8)
ax.contour(X_grid, Y_grid, p_x_given_1, levels=5, colors='C1', alpha=0.8)
ax.scatter(*mu0, s=100, c='C0', marker='x', linewidths=3, zorder=5)
ax.scatter(*mu1, s=100, c='C1', marker='x', linewidths=3, zorder=5)
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('Distribution vraie $p(x, y)$\n(inconnue)')
ax.set_aspect('equal')

# Right: Finite sample (what we observe)
ax = axes[1]
ax.scatter(X[y == 0, 0], X[y == 0, 1], c='C0', alpha=0.7, s=50, label='Classe 0')
ax.scatter(X[y == 1, 0], X[y == 1, 1], c='C1', alpha=0.7, s=50, label='Classe 1')
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title(f'Échantillon observé $\\mathcal{{D}}$\n($N = {len(X)}$ points)')
ax.legend()
ax.set_aspect('equal')

plt.tight_layout()
<Figure size 1200x500 with 2 Axes>

Nous ne voyons que les points à droite. La structure continue à gauche, incluant les contours, les densités, ainsi que les régions de haute et basse probabilité, nous est cachée. C’est à partir de ces quelques points que nous devons estimer la performance de notre modèle.

Obstacle 2: L’intégration est computationnellement intractable

Supposons, par un miracle, que nous connaissions exactement p(x,y)p(\mathbf{x}, y). Pourrions-nous alors calculer le risque R(f)=(y,f(x))p(x,y)dxdy\mathcal{R}(f) = \int \ell(y, f(\mathbf{x})) \, p(\mathbf{x}, y) \, d\mathbf{x} \, dy?

La réponse est généralement non, pour plusieurs raisons:

Pour les espaces continus: L’intégrale est une intégrale de grande dimension. Si xRd\mathbf{x} \in \mathbb{R}^d avec dd grand (par exemple, d=1000d = 1000 pour des images ou d=106d = 10^6 pour des données textuelles), nous devons intégrer sur un espace de dimension d+1d+1.

Pour vous rappeler l’idée, en calcul on approche une intégrale en 1D par une somme: on découpe l’intervalle en petites tranches et on additionne des aires de rectangles ou de trapèzes. Par exemple, sur [a,b][a,b]:

abg(x)dx    m=1Mg(xm)Δx\int_a^b g(x)\,dx \;\approx\; \sum_{m=1}^{M} g(x_m)\,\Delta x

Cette idée générale, qui consiste à remplacer une intégrale par une somme pondérée de valeurs de gg évaluées à des points xmx_m, s’appelle l’intégration numérique (ou quadrature).

Le problème en apprentissage est que notre intégrale n’est pas en 1D. Si on applique le même raisonnement en dimension dd en mettant, disons, MM points par dimension, on obtient une grille de taille MdM^d (et ici dd peut être très grand). Le nombre de points à évaluer explose donc exponentiellement avec dd. C’est exactement la malédiction de la dimensionnalité.

Pour les espaces discrets: Si x\mathbf{x} et yy sont discrets mais prennent de nombreuses valeurs, la somme xy(y,f(x))p(x,y)\sum_{\mathbf{x}} \sum_y \ell(y, f(\mathbf{x})) \cdot p(\mathbf{x}, y) peut avoir un nombre exponentiel de termes. Par exemple, si x\mathbf{x} est un vecteur binaire de dimension dd, il y a 2d2^d valeurs possibles pour x\mathbf{x}. Pour d=100d = 100, cela fait déjà 210010302^{100} \approx 10^{30} termes à sommer, ce qui est computationnellement impossible.

Intégration de Monte Carlo: On pourrait penser utiliser l’intégration de Monte Carlo: tirer des échantillons (x,y)(\mathbf{x}, y) selon p(x,y)p(\mathbf{x}, y) et estimer l’intégrale par la moyenne empirique. Mais pour obtenir une estimation précise du risque, nous aurions besoin d’un très grand nombre d’échantillons (potentiellement infini pour une précision parfaite). De plus, cela nécessiterait de pouvoir échantillonner efficacement depuis p(x,y)p(\mathbf{x}, y), ce qui est lui-même un problème difficile si la distribution est complexe.

La figure suivante illustre la malédiction de la dimensionnalité. Avec seulement 10 points par dimension pour une quadrature numérique, le nombre total de points d’évaluation explose rapidement.

Source
import numpy as np
import matplotlib.pyplot as plt

# Number of grid points per dimension
points_per_dim = 10

# Dimensions to consider
dimensions = np.array([1, 2, 3, 5, 10, 20, 50, 100])

# Total grid points = points_per_dim^d
total_points = points_per_dim ** dimensions.astype(float)

fig, ax = plt.subplots(figsize=(8, 5))

bars = ax.bar(range(len(dimensions)), total_points, color='steelblue', edgecolor='black')

# Add reference lines
ax.axhline(y=1e9, color='C1', linestyle='--', alpha=0.7, label='1 milliard (limite pratique)')
ax.axhline(y=1e80, color='C3', linestyle=':', alpha=0.7, label='$10^{80}$ (atomes dans l\'univers)')

ax.set_yscale('log')
ax.set_xticks(range(len(dimensions)))
ax.set_xticklabels([f'd={d}' for d in dimensions])
ax.set_xlabel('Dimension de l\'espace des entrées')
ax.set_ylabel('Nombre de points de grille')
ax.set_title(f'Points nécessaires pour l\'intégration numérique\n({points_per_dim} points par dimension)')
ax.legend(loc='upper left')

# Annotate a few bars
for i, (d, n) in enumerate(zip(dimensions, total_points)):
    if d <= 5:
        ax.annotate(f'$10^{{{d}}}$', (i, n), ha='center', va='bottom', fontsize=9)
    elif d == 10:
        ax.annotate(f'$10^{{{10}}}$', (i, n), ha='center', va='bottom', fontsize=9)
    elif d == 100:
        ax.annotate(f'$10^{{{100}}}$', (i, n), ha='center', va='bottom', fontsize=9)

ax.set_ylim(1, 1e105)

plt.tight_layout()
<Figure size 800x500 with 1 Axes>

En dimension 10, il faut déjà 1010 points, soit dix milliards. En dimension 100, il en faut 10100, un nombre qui dépasse le nombre d’atomes dans l’univers observable. L’intégration numérique directe est donc impossible pour les problèmes de haute dimension, même si nous connaissions p(x,y)p(\mathbf{x}, y) exactement.

Le risque empirique comme seule option pratique

Face à ces obstacles, le risque empirique est notre seule option calculable. Mais il y a une bonne nouvelle: le risque empirique est une forme d’intégration de Monte Carlo, et Monte Carlo a une propriété remarquable.

MéthodeComplexitéExigence
Quadrature (règles trapézoïdales, etc.)O(Md)O(M^d)Connaître p(x,y)p(\mathbf{x},y) exactement
Monte CarloO(N)O(N)Avoir des échantillons de p(x,y)p(\mathbf{x},y)

La complexité de Monte Carlo est indépendante de la dimension dd. Elle ne dépend que du nombre d’échantillons NN. C’est cette propriété qui rend l’apprentissage possible en haute dimension. De plus, nous n’avons pas besoin de connaître la valeur numérique de p(x,y)p(\mathbf{x},y): nous avons seulement besoin de pouvoir tirer des échantillons de cette distribution. C’est exactement ce que nos données d’entraînement nous fournissent.

Le risque empirique remplace l’intégrale sur la distribution inconnue par une moyenne sur l’échantillon fini que nous possédons:

R^(f,D)=1Ni=1N(yi,f(xi))\hat{\mathcal{R}}(f, \mathcal{D}) = \frac{1}{N} \sum_{i=1}^{N} \ell(y_i, f(\mathbf{x}_i))

Cette formule est directe à évaluer: nous parcourons nos NN exemples d’entraînement, calculons la perte pour chacun, et faisons la moyenne.

Reprenons les données de freinage. Divisons-les en deux parties: les mesures à vitesses faibles (4-19 mph) pour l’entraînement, et les mesures à vitesses élevées (20-25 mph) pour le test. Le risque empirique sur l’ensemble d’entraînement mesure la qualité de l’ajustement. Le risque empirique sur l’ensemble de test estime la performance sur des vitesses non vues pendant l’entraînement.

Source
import numpy as np
import matplotlib.pyplot as plt

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

# Split: train on low speeds, test on high speeds
train_mask = speed < 20
test_mask = speed >= 20

speed_train, dist_train = speed[train_mask], dist[train_mask]
speed_test, dist_test = speed[test_mask], dist[test_mask]

# Fit on training data
coeffs = np.polyfit(speed_train, dist_train, 2)

# Compute MSE on train and test
pred_train = np.polyval(coeffs, speed_train)
pred_test = np.polyval(coeffs, speed_test)
mse_train = np.mean((dist_train - pred_train)**2)
mse_test = np.mean((dist_test - pred_test)**2)

plt.figure(figsize=(7, 4))
plt.scatter(speed_train, dist_train, alpha=0.7, label=f'Entraînement (EQM={mse_train:.1f})')
plt.scatter(speed_test, dist_test, alpha=0.7, marker='s', label=f'Test (EQM={mse_test:.1f})')

speed_grid = np.linspace(4, 28, 100)
plt.plot(speed_grid, np.polyval(coeffs, speed_grid), 'k--', alpha=0.6, label='Fonction ajustée')

plt.axvline(x=20, color='gray', linestyle=':', alpha=0.5)
plt.xlabel('Vitesse (mph)')
plt.ylabel('Distance (ft)')
plt.legend()
plt.tight_layout()
<Figure size 700x400 with 1 Axes>

Dans cet exemple, l’EQM (MSE, erreur quadratique moyenne) sur l’ensemble de test est plus élevé que sur l’ensemble d’entraînement. Cet écart est typique: la fonction a été optimisée pour les données d’entraînement, pas pour les données de test.

Sous l’hypothèse que les exemples (xi,yi)(\mathbf{x}_i, y_i) sont tirés indépendamment et identiquement distribués (i.i.d.) selon p(x,y)p(\mathbf{x}, y), le risque empirique est un estimateur non biaisé du vrai risque: E[R^(f,D)]=R(f)\mathbb{E}[\hat{\mathcal{R}}(f, \mathcal{D})] = \mathcal{R}(f). Cela signifie qu’en moyenne, sur tous les échantillons possibles, le risque empirique est égal au vrai risque.

Par la loi des grands nombres, lorsque NN \to \infty, le risque empirique converge vers le vrai risque (presque sûrement). Avec suffisamment de données, si l’échantillon est représentatif de la distribution, le risque empirique devrait être proche du risque.

La figure suivante illustre cette convergence. Nous utilisons le problème de classification gaussienne pour lequel nous pouvons calculer le vrai risque analytiquement. Chaque courbe montre l’évolution du risque empirique pour un échantillon de taille croissante. Toutes les courbes convergent vers le vrai risque (ligne pointillée), mais avec des fluctuations qui diminuent à mesure que NN augmente.

Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Paramètres du mélange gaussien
mu0, mu1 = np.array([0.0, 0.0]), np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])

def gaussian_pdf(x, mu, cov):
    d = len(mu)
    diff = x - mu
    cov_inv = np.linalg.inv(cov)
    mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
    return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))

# Compute Bayes-optimal classifier error rate (true risk)
# For Gaussian classes with equal covariance, the Bayes error is:
# P(error) = Phi(-d/2) where d is the Mahalanobis distance between means
cov_inv = np.linalg.inv(cov)
d_squared = (mu1 - mu0) @ cov_inv @ (mu1 - mu0)
d = np.sqrt(d_squared)
true_risk = norm.cdf(-d / 2)

# Simulate empirical risk for different sample sizes
sample_sizes = np.arange(10, 1001, 10)
n_runs = 20

fig, ax = plt.subplots(figsize=(9, 5))

# Store all runs for confidence band
all_risks = np.zeros((n_runs, len(sample_sizes)))

for run in range(n_runs):
    empirical_risks = []
    # Generate a large dataset and compute cumulative empirical risk
    rng = np.random.default_rng(run)
    n_total = 1000
    X0 = rng.multivariate_normal(mu0, cov, n_total // 2)
    X1 = rng.multivariate_normal(mu1, cov, n_total // 2)
    X = np.vstack([X0, X1])
    y = np.concatenate([np.zeros(n_total // 2), np.ones(n_total // 2)])
    perm = rng.permutation(n_total)
    X, y = X[perm], y[perm]
    
    # Bayes-optimal classifier: predict 1 if w^T x + b > 0
    w = cov_inv @ (mu1 - mu0)
    b = -0.5 * (mu1 @ cov_inv @ mu1 - mu0 @ cov_inv @ mu0)
    
    for n in sample_sizes:
        X_n, y_n = X[:n], y[:n]
        predictions = (X_n @ w + b > 0).astype(float)
        emp_risk = np.mean(predictions != y_n)
        empirical_risks.append(emp_risk)
    
    all_risks[run] = empirical_risks
    ax.plot(sample_sizes, empirical_risks, 'C0-', alpha=0.15, linewidth=0.8)

# Mean and confidence bands
mean_risk = np.mean(all_risks, axis=0)
std_risk = np.std(all_risks, axis=0)
ax.fill_between(sample_sizes, mean_risk - 2*std_risk, mean_risk + 2*std_risk, 
                alpha=0.3, color='C0', label='Intervalle ± 2 écarts-types')
ax.plot(sample_sizes, mean_risk, 'C0-', linewidth=2, label='Moyenne empirique')

# True risk
ax.axhline(y=true_risk, color='C3', linestyle='--', linewidth=2, 
           label=f'Vrai risque = {true_risk:.3f}')

ax.set_xlabel('Taille de l\'échantillon $N$')
ax.set_ylabel('Risque empirique (taux d\'erreur)')
ax.set_title('Convergence du risque empirique vers le vrai risque')
ax.legend(loc='upper right')
ax.set_xlim(0, 1000)
ax.set_ylim(0, 0.35)

plt.tight_layout()
<Figure size 900x500 with 1 Axes>

Avec N=50N = 50, le risque empirique peut facilement varier de 0.10 à 0.25 selon l’échantillon. Avec N=500N = 500, la variabilité est beaucoup plus faible. C’est la loi des grands nombres en action: plus l’échantillon est grand, plus l’estimation est précise.

Le compromis fondamental

Cette situation crée un compromis fondamental en apprentissage automatique:

L’écart entre ces deux quantités est au cœur de l’apprentissage automatique. Un modèle peut avoir un risque empirique très faible (il performe bien sur les données d’entraînement) tout en ayant un risque élevé (il performe mal sur de nouvelles données). C’est le problème du surapprentissage, que nous explorerons plus en détail dans le chapitre sur la généralisation.

La question de savoir quand et à quelle vitesse l’approximation du risque par le risque empirique est fiable relève de la théorie de la généralisation, que nous aborderons au chapitre suivant.

Minimisation du risque empirique

Nous avons maintenant les éléments pour formuler l’apprentissage comme un problème d’optimisation. Nous cherchons la fonction ff dans une classe F\mathcal{F} qui minimise le risque:

f=argminfFR(f)f^\star = \arg\min_{f \in \mathcal{F}} \mathcal{R}(f)

Puisque le risque est inaccessible, nous le remplaçons par le risque empirique:

f^=argminfFR^(f,D)\hat{f} = \arg\min_{f \in \mathcal{F}} \hat{\mathcal{R}}(f, \mathcal{D})

Ce principe est la minimisation du risque empirique (MRE): choisir la fonction qui fait le moins d’erreurs sur les données d’entraînement, en espérant que cette performance se transfère aux nouvelles données.

La classe F\mathcal{F} est notre classe d’hypothèses. Elle représente l’ensemble des fonctions que nous sommes prêts à considérer. Le choix de F\mathcal{F} encode nos hypothèses sur la forme de la relation entre entrées et sorties.

Un premier exemple: les modèles linéaires

Pour rendre ces concepts concrets, considérons la classe la plus simple: les modèles linéaires. Un modèle linéaire suppose que la sortie est une combinaison linéaire des entrées:

f(x;θ)=θ0+j=1dθjxj=θxf(\mathbf{x}; \boldsymbol{\theta}) = \theta_0 + \sum_{j=1}^d \theta_j x_j = \boldsymbol{\theta}^\top \mathbf{x}

xRd+1\mathbf{x} \in \mathbb{R}^{d+1} est le vecteur d’entrée augmenté d’un 1 pour le biais (x0=1x_0 = 1), et θRd+1\boldsymbol{\theta} \in \mathbb{R}^{d+1} est le vecteur de paramètres contenant le biais θ0\theta_0 et les poids θ1,,θd\theta_1, \ldots, \theta_d.

Cette forme est restrictive: elle suppose que la relation entre entrées et sorties est linéaire. Pour les données de freinage, cela signifierait que la distance est proportionnelle à la vitesse, ce qui n’est pas le cas (la relation est plutôt quadratique). Néanmoins, les modèles linéaires sont utiles comme point de départ, et nous verrons comment les étendre pour capturer des relations non linéaires.

Avec cette classe F\mathcal{F} fixée, l’apprentissage consiste à trouver les paramètres θ\boldsymbol{\theta} qui minimisent le risque empirique. Pour la perte quadratique, cela revient à minimiser la somme des carrés des résidus.

Mais quand le minimiseur du risque empirique a-t-il un faible risque? Si f^\hat{f} minimise R^\hat{\mathcal{R}} et ff^\star minimise R\mathcal{R}, nous voulons que R(f^)\mathcal{R}(\hat{f}) soit proche de R(f)\mathcal{R}(f^\star). La réponse dépend de la taille de l’échantillon NN, de la complexité de la classe F\mathcal{F}, et de propriétés de la distribution pp.

Résoudre le problème d’optimisation

Quand nous écrivons argminθ\arg\min_{\boldsymbol{\theta}}, nous cherchons les paramètres qui rendent la fonction objectif aussi petite que possible. Mais comment trouver ces paramètres en pratique?

La réponse dépend de la forme du problème:

Exemple: solution analytique pour la régression linéaire (MCO)

Pour illustrer les solutions analytiques, considérons la régression linéaire avec perte quadratique. L’objectif est de minimiser la somme des carrés des résidus:

RSS(θ)=i=1N(yiθxi)2=yXθ22\text{RSS}(\boldsymbol{\theta}) = \sum_{i=1}^N (y_i - \boldsymbol{\theta}^\top \mathbf{x}_i)^2 = \|\mathbf{y} - \mathbf{X}\boldsymbol{\theta}\|_2^2

X\mathbf{X} est la matrice N×(d+1)N \times (d+1) des entrées (avec une colonne de 1 pour le biais) et y\mathbf{y} est le vecteur des sorties.

En développant et en calculant le gradient:

θRSS(θ)=2Xy+2XXθ\nabla_{\boldsymbol{\theta}} \text{RSS}(\boldsymbol{\theta}) = -2\mathbf{X}^\top \mathbf{y} + 2\mathbf{X}^\top \mathbf{X} \boldsymbol{\theta}

En posant le gradient égal à zéro, nous obtenons les équations normales:

XXθ=Xy\mathbf{X}^\top \mathbf{X} \boldsymbol{\theta} = \mathbf{X}^\top \mathbf{y}

Si la matrice XX\mathbf{X}^\top \mathbf{X} est inversible, la solution unique est:

θ^MCO=(XX)1Xy\hat{\boldsymbol{\theta}}_{\text{MCO}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}

Cette solution porte le nom d’estimateur des moindres carrés ordinaires (MCO, ou ordinary least squares, OLS). Elle peut être calculée directement sans itération, ce qui en fait un exemple classique de solution analytique.

Le chapitre sur l’optimisation présentera ces méthodes en détail. Pour l’instant, nous nous concentrons sur la formulation des problèmes d’apprentissage comme problèmes d’optimisation, en gardant à l’esprit que des outils existent pour les résoudre.

Généralisation

La différence entre le risque et le risque empirique est l’écart de généralisation:

Eˊcart=R(f)R^(f;Dtrain)\text{Écart} = \mathcal{R}(f) - \hat{\mathcal{R}}(f; \mathcal{D}_{\text{train}})

Un modèle qui minimise le risque empirique peut avoir un risque élevé si cet écart est grand. Ce phénomène est le surapprentissage: le modèle s’ajuste aux particularités de l’échantillon d’entraînement, y compris le bruit, plutôt qu’aux régularités sous-jacentes. L’erreur d’entraînement est faible, mais l’erreur sur de nouvelles données est élevée.

À l’inverse, un modèle trop simple peut avoir un risque empirique et un risque tous deux élevés. C’est le sous-apprentissage: le modèle n’a pas la capacité de capturer la structure des données.

Extrapolation

Un cas particulier de mauvaise généralisation est l’extrapolation: prédire pour des entrées en dehors de la plage des données d’entraînement. Même un modèle bien ajusté peut échouer spectaculairement lorsqu’on lui demande de prédire au-delà de ce qu’il a vu.

Considérons des essais en soufflerie pour mesurer la portance d’une aile à différentes vitesses. Les tests sont effectués entre 20 et 60 m/s. L’ingénieur veut prédire la portance à 100 m/s.

Source
import numpy as np
import matplotlib.pyplot as plt

# Données de portance aérodynamique (simulées)
np.random.seed(42)
rho, S, C_L = 1.225, 20.0, 0.5
v_train = np.linspace(20, 60, 8)
L_true_train = 0.5 * rho * v_train**2 * S * C_L
L_train = L_true_train + np.random.normal(0, 400, len(v_train))

coeffs_2 = np.polyfit(v_train, L_train, 2)
coeffs_5 = np.polyfit(v_train, L_train, 5)

v_extrap = np.linspace(15, 110, 200)
L_true_extrap = 0.5 * rho * v_extrap**2 * S * C_L

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

for ax, coeffs, deg in zip(axes, [coeffs_2, coeffs_5], [2, 5]):
    ax.scatter(v_train, L_train, s=50, zorder=5, label='Observations')
    ax.plot(v_extrap, L_true_extrap, 'g-', alpha=0.5, label='Vraie relation')
    L_pred = np.polyval(coeffs, v_extrap)
    ax.plot(v_extrap, L_pred, 'k--', label=f'Polynôme degré {deg}')
    ax.axvline(60, color='gray', linestyle=':', alpha=0.5)
    ax.axvspan(60, 110, alpha=0.1, color='red')
    ax.set_xlabel('Vitesse (m/s)')
    ax.set_ylabel('Portance (N)')
    ax.set_title(f'Degré {deg}')
    ax.legend(loc='upper left')
    ax.set_ylim(-5000, 80000)
    ax.text(85, 5000, 'Extrapolation', ha='center', fontsize=10, color='red', alpha=0.7)

plt.tight_layout()
<Figure size 1000x400 with 2 Axes>

Le polynôme de degré 2 (qui correspond au vrai modèle physique Lv2L \propto v^2) extrapole correctement. Le polynôme de degré 5, bien qu’il ajuste aussi bien les données d’entraînement, diverge complètement en dehors de la plage observée.

Régularisation

Une manière de contrôler le surapprentissage consiste à pénaliser la complexité du modèle directement dans la fonction objectif. Le risque empirique régularisé est:

R^λ(θ)=R^(θ)+λC(θ)\hat{\mathcal{R}}_\lambda(\boldsymbol{\theta}) = \hat{\mathcal{R}}(\boldsymbol{\theta}) + \lambda \, C(\boldsymbol{\theta})

C(θ)C(\boldsymbol{\theta}) mesure la complexité du modèle et λ0\lambda \geq 0 contrôle l’intensité de la pénalisation. Un choix courant est la régularisation 2\ell_2 (ou weight decay):

C(θ)=θ22=jθj2C(\boldsymbol{\theta}) = \|\boldsymbol{\theta}\|_2^2 = \sum_j \theta_j^2

Cette pénalisation pousse les paramètres vers zéro, ce qui a pour effet de lisser la fonction apprise. En régression linéaire, l’ajout de cette pénalité donne la régression ridge:

θ^ridge=argminθ1Ni=1N(yiθxi)2+λθ22\hat{\boldsymbol{\theta}}_{\text{ridge}} = \arg\min_{\boldsymbol{\theta}} \frac{1}{N}\sum_{i=1}^N (y_i - \boldsymbol{\theta}^\top \mathbf{x}_i)^2 + \lambda \|\boldsymbol{\theta}\|_2^2

Illustrons l’effet de la régularisation sur le même problème de régression polynomiale. Avec un polynôme de degré 15 et différentes valeurs de λ\lambda:

Source
import numpy as np
import matplotlib.pyplot as plt

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

# Train/test split
np.random.seed(42)
indices = np.random.permutation(len(speed))
train_idx, test_idx = indices[:35], indices[35:]
speed_train, dist_train = speed[train_idx], dist[train_idx]
speed_test, dist_test = speed[test_idx], dist[test_idx]

# Build polynomial features (degree 15)
degree = 15
def poly_features(x, deg):
    return np.vstack([x**i for i in range(deg+1)]).T

X_train = poly_features(speed_train, degree)
X_test = poly_features(speed_test, degree)

# Ridge regression for different lambda values
lambdas = [0, 1e-6, 1e-3, 1]
fig, axes = plt.subplots(2, 2, figsize=(10, 8))

for ax, lam in zip(axes.flat, lambdas):
    # Ridge solution: (X^T X + lambda I)^{-1} X^T y
    I = np.eye(X_train.shape[1])
    I[0, 0] = 0  # Don't regularize bias
    w = np.linalg.solve(X_train.T @ X_train + lam * I, X_train.T @ dist_train)
    
    # Predictions
    pred_train = X_train @ w
    pred_test = X_test @ w
    mse_train = np.mean((dist_train - pred_train)**2)
    mse_test = np.mean((dist_test - pred_test)**2)
    
    # Plot
    ax.scatter(speed_train, dist_train, alpha=0.6, s=30, label='Entraînement')
    ax.scatter(speed_test, dist_test, alpha=0.6, s=30, marker='s', label='Test')
    
    speed_grid = np.linspace(3, 26, 200)
    X_grid = poly_features(speed_grid, degree)
    pred_grid = X_grid @ w
    pred_grid = np.clip(pred_grid, -50, 200)
    ax.plot(speed_grid, pred_grid, 'k-', alpha=0.7)
    
    ax.set_xlim(3, 26)
    ax.set_ylim(-20, 150)
    ax.set_xlabel('Vitesse (mph)')
    ax.set_ylabel('Distance (ft)')
    ax.set_title(f'$\\lambda$ = {lam}: Entr. EQM={mse_train:.1f}, Test EQM={mse_test:.1f}')
    if lam == 0:
        ax.legend()

plt.tight_layout()
<Figure size 1000x800 with 4 Axes>

Sans régularisation (λ=0\lambda = 0), le polynôme de degré 15 oscille fortement. Avec une régularisation modérée (λ=103\lambda = 10^{-3}), les oscillations sont atténuées et l’erreur de test diminue. Avec une régularisation trop forte (λ=1\lambda = 1), le modèle devient trop contraint et sous-apprend.

Solution analytique de la régression ridge

Comme pour les moindres carrés ordinaires, la régression ridge admet une solution analytique. L’objectif régularisé est:

RSSλ(θ)=yXθ22+λθ22\text{RSS}_\lambda(\boldsymbol{\theta}) = \|\mathbf{y} - \mathbf{X}\boldsymbol{\theta}\|_2^2 + \lambda \|\boldsymbol{\theta}\|_2^2

En développant et en calculant le gradient:

θRSSλ(θ)=2Xy+2XXθ+2λθ=2Xy+2(XX+λI)θ\nabla_{\boldsymbol{\theta}} \text{RSS}_\lambda(\boldsymbol{\theta}) = -2\mathbf{X}^\top \mathbf{y} + 2\mathbf{X}^\top \mathbf{X} \boldsymbol{\theta} + 2\lambda \boldsymbol{\theta} = -2\mathbf{X}^\top \mathbf{y} + 2(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}) \boldsymbol{\theta}

En posant le gradient égal à zéro, nous obtenons les équations normales régularisées:

(XX+λI)θ=Xy(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}) \boldsymbol{\theta} = \mathbf{X}^\top \mathbf{y}

La solution est:

θ^ridge=(XX+λI)1Xy\hat{\boldsymbol{\theta}}_{\text{ridge}} = (\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X}^\top \mathbf{y}

Comparons avec la solution MCO: θ^MCO=(XX)1Xy\hat{\boldsymbol{\theta}}_{\text{MCO}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}. La seule différence est l’ajout du terme λI\lambda \mathbf{I} à la matrice XX\mathbf{X}^\top \mathbf{X}.

Solution via décomposition en valeurs singulières (SVD)

Cette section présente une autre façon d’exprimer les solutions MCO et Ridge, en utilisant la décomposition en valeurs singulières (SVD). Si vous n’avez jamais rencontré la SVD, ne vous inquiétez pas: nous allons l’introduire progressivement. Cette approche n’est pas strictement nécessaire pour comprendre la régression ridge, mais elle offre une interprétation géométrique très éclairante qui révèle pourquoi la régularisation fonctionne.

Qu’est-ce que la SVD? Si vous avez déjà rencontré la décomposition en valeurs propres, la SVD en est une généralisation. Pour une matrice carrée symétrique A\mathbf{A}, la décomposition en valeurs propres s’écrit A=QΛQ\mathbf{A} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^\top, où Q\mathbf{Q} contient les vecteurs propres et Λ\boldsymbol{\Lambda} les valeurs propres. La SVD généralise cette idée à n’importe quelle matrice, même rectangulaire.

Pour une matrice X\mathbf{X} de données, la SVD la réécrit comme le produit de trois matrices:

X=UDV\mathbf{X} = \mathbf{U} \mathbf{D} \mathbf{V}^\top

Lien avec la décomposition en valeurs propres: Les colonnes de V\mathbf{V} sont les vecteurs propres de XX\mathbf{X}^\top \mathbf{X}, et les valeurs singulières djd_j sont les racines carrées des valeurs propres de XX\mathbf{X}^\top \mathbf{X}. Autrement dit, si XX=VΛV\mathbf{X}^\top \mathbf{X} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^\top est la décomposition en valeurs propres de XX\mathbf{X}^\top \mathbf{X}, alors dj=λjd_j = \sqrt{\lambda_j}λj\lambda_j sont les valeurs propres. De même, les colonnes de U\mathbf{U} sont les vecteurs propres de XX\mathbf{X} \mathbf{X}^\top.

Cette connexion est utile car XX\mathbf{X}^\top \mathbf{X} apparaît naturellement dans la régression (c’est la matrice que nous inversons pour MCO). Les valeurs singulières djd_j nous renseignent donc directement sur le “conditionnement” de cette matrice: si certaines valeurs singulières sont très petites, alors XX\mathbf{X}^\top \mathbf{X} est proche d’être singulière (non inversible).

Interprétation géométrique:

Solution MCO via SVD: En utilisant cette décomposition, la solution MCO peut s’écrire:

θ^MCO=j=1dujydjvj\hat{\boldsymbol{\theta}}_{\text{MCO}} = \sum_{j=1}^d \frac{\mathbf{u}_j^\top \mathbf{y}}{d_j} \mathbf{v}_j

Cette formule décompose la solution en une somme de contributions le long de chaque direction principale vj\mathbf{v}_j. Le terme ujydj\frac{\mathbf{u}_j^\top \mathbf{y}}{d_j} mesure combien la sortie y\mathbf{y} s’aligne avec la direction uj\mathbf{u}_j, divisé par l’amplitude djd_j de cette direction. Notez que diviser par une petite valeur singulière djd_j peut amplifier le bruit, ce qui explique pourquoi MCO peut être instable quand certaines directions ont de petites valeurs singulières.

Solution Ridge via SVD: Pour Ridge, la solution devient:

θ^ridge=j=1ddj2dj2+λujydjvj\hat{\boldsymbol{\theta}}_{\text{ridge}} = \sum_{j=1}^d \frac{d_j^2}{d_j^2 + \lambda} \frac{\mathbf{u}_j^\top \mathbf{y}}{d_j} \mathbf{v}_j

La différence avec MCO est le facteur de rétrécissement dj2dj2+λ\frac{d_j^2}{d_j^2 + \lambda} qui multiplie chaque terme. Ce facteur est toujours inférieur à 1, ce qui “rétrécit” chaque composante vers zéro. L’effet clé est que ce rétrécissement est différencié: les directions avec de petites valeurs singulières sont rétrécies plus fortement que celles avec de grandes valeurs singulières.

Interprétation géométrique: Les directions principales vj\mathbf{v}_j définissent les axes d’une ellipse de confiance dans l’espace des coefficients. Pour MCO, les longueurs des demi-axes sont proportionnelles à 1/dj1/d_j: plus une direction a une petite valeur singulière djd_j, plus la variance d’estimation est grande (nous verrons cette distinction cruciale plus loin). Avec Ridge, ces longueurs deviennent proportionnelles à djdj2+λ\frac{d_j}{d_j^2 + \lambda}: l’ellipse se rétrécit globalement, et plus rapidement le long des directions associées aux petites valeurs singulières. C’est exactement ce que nous voulons: réduire la variance d’estimation là où elle est la plus grande, c’est-à-dire dans les directions où les données sont peu dispersées.

Avantages numériques: Au-delà de l’interprétation, la SVD offre aussi des avantages pratiques. Elle est plus stable numériquement que l’inversion directe de XX\mathbf{X}^\top \mathbf{X}, surtout quand cette matrice est mal conditionnée (c’est-à-dire quand certaines valeurs singulières sont très petites). Les algorithmes SVD gèrent mieux ces cas délicats.

Visualisation: ellipse des données et vecteurs singuliers

Pour rendre ces concepts concrets, visualisons ce que la SVD capture sur un nuage de données 2D. Générons des points suivant une distribution gaussienne avec une covariance non triviale (les deux variables sont corrélées).

Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

np.random.seed(42)

# Générer des données gaussiennes corrélées
n_points = 200
mean = [0, 0]
cov = [[2.0, 1.2], [1.2, 1.0]]  # Covariance non diagonale
X = np.random.multivariate_normal(mean, cov, n_points)

# Centrer les données
X_centered = X - X.mean(axis=0)

# SVD de la matrice de données centrée
U, d, Vt = np.linalg.svd(X_centered, full_matrices=False)
V = Vt.T

# Les valeurs singulières sont liées aux écarts-types: d_j / sqrt(N-1)
# Pour l'ellipse, nous utilisons les écarts-types dans chaque direction
std_1 = d[0] / np.sqrt(n_points - 1)
std_2 = d[1] / np.sqrt(n_points - 1)

# Créer la figure
fig, ax = plt.subplots(figsize=(8, 6))

# Tracer les points
ax.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.5, s=20, c='tab:blue', label='Données')

# Tracer les vecteurs singuliers (directions principales)
origin = [0, 0]
scale = 2  # Facteur d'échelle pour la visualisation

# Premier vecteur singulier (direction de plus grande variance)
ax.annotate('', xy=V[:, 0] * std_1 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:red', lw=2.5))
ax.annotate('', xy=-V[:, 0] * std_1 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:red', lw=2.5))

# Deuxième vecteur singulier (direction de plus petite variance)
ax.annotate('', xy=V[:, 1] * std_2 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:orange', lw=2.5))
ax.annotate('', xy=-V[:, 1] * std_2 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:orange', lw=2.5))

# Ellipse de confiance (2 écarts-types)
angle = np.degrees(np.arctan2(V[1, 0], V[0, 0]))
ellipse = Ellipse(xy=(0, 0), width=4*std_1, height=4*std_2, angle=angle,
                  fill=False, edgecolor='gray', linestyle='--', linewidth=1.5)
ax.add_patch(ellipse)

# Annotations
ax.text(V[0, 0] * std_1 * scale * 1.15, V[1, 0] * std_1 * scale * 1.15, 
        f'$\\mathbf{{v}}_1$ ($d_1 = {d[0]:.1f}$)', fontsize=11, color='tab:red')
ax.text(V[0, 1] * std_2 * scale * 1.3, V[1, 1] * std_2 * scale * 1.3, 
        f'$\\mathbf{{v}}_2$ ($d_2 = {d[1]:.1f}$)', fontsize=11, color='tab:orange')

ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('Nuage gaussien et directions principales (SVD)')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)
ax.set_xlim(-4, 4)
ax.set_ylim(-3, 3)

plt.tight_layout()
<Figure size 800x600 with 1 Axes>

La figure montre un nuage de 200 points tirés d’une gaussienne 2D. Les flèches représentent les vecteurs singuliers v1\mathbf{v}_1 et v2\mathbf{v}_2:

L’ellipse en pointillés représente la région contenant environ 95% des données si elles suivent exactement la distribution gaussienne. Ses axes coïncident avec les vecteurs singuliers, et les longueurs des demi-axes sont proportionnelles aux valeurs singulières.

Cette visualisation illustre pourquoi la SVD est si utile: elle identifie automatiquement les axes naturels des données. Dans le contexte de la régression, si les caractéristiques forment un nuage allongé (valeurs singulières très différentes), alors certaines directions contiennent beaucoup d’information (grandes valeurs singulières) tandis que d’autres en contiennent peu (petites valeurs singulières).

Deux variances: données vs estimation

Avant d’aller plus loin, clarifions une source fréquente de confusion. Le mot variance apparaît dans deux contextes très différents lorsqu’on parle de SVD et de régularisation:

  1. Variance des données (dispersion): mesure l’étalement des données le long d’une direction vj\mathbf{v}_j. Elle est proportionnelle à dj2d_j^2. Une grande valeur singulière djd_j signifie que les données sont très dispersées dans cette direction.

  2. Variance d’estimation (incertitude): mesure l’incertitude sur notre estimé θ^j\hat{\theta}_j du paramètre correspondant à la direction vj\mathbf{v}_j. Elle est proportionnelle à 1/dj21/d_j^2. Une petite valeur singulière djd_j signifie que notre estimé est très incertain.

Ces deux variances sont inversement reliées:

Valeur singulière djd_jVariance des donnéesVariance d’estimationInterprétation
GrandeÉlevée (données étalées)Faible (estimé précis)Beaucoup d’information
PetiteFaible (données concentrées)Élevée (estimé incertain)Peu d’information

Intuition: Imaginez estimer une pente à partir de données. Si les points sont très étalés horizontalement (grande variance des données en xx), la pente est facile à déterminer avec précision (faible variance d’estimation). Si les points sont tous regroupés (petite variance des données), la pente est très incertaine (grande variance d’estimation).

C’est cette relation inverse qui explique le comportement de Ridge:

Ainsi, quand nous disons que « Ridge contrôle la variance », nous parlons de la variance d’estimation des paramètres, pas de la variance des données. La régularisation n’affecte pas la dispersion des données; elle réduit l’incertitude de nos estimés en les « tirant » vers zéro.

Spectre des valeurs singulières et rang effectif

En pratique, les données réelles ont souvent des dizaines ou des centaines de dimensions. Comment se comportent les valeurs singulières dans ce cas? Examinons un exemple avec des données de dimension plus élevée.

Source
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Simuler des données avec une structure de rang bas + bruit
n_samples = 100
n_features = 30

# Vraie structure: combinaison de 5 facteurs latents
n_latent = 5
latent = np.random.randn(n_samples, n_latent)
loadings = np.random.randn(n_latent, n_features)
X_signal = latent @ loadings

# Ajouter du bruit
noise_level = 0.5
X = X_signal + noise_level * np.random.randn(n_samples, n_features)

# Centrer
X_centered = X - X.mean(axis=0)

# SVD
U, d, Vt = np.linalg.svd(X_centered, full_matrices=False)

# Variance expliquée
variance_explained = d**2 / np.sum(d**2)
cumulative_variance = np.cumsum(variance_explained)

# Créer la figure
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Panneau gauche: spectre des valeurs singulières
ax1 = axes[0]
ax1.semilogy(range(1, len(d)+1), d, 'o-', markersize=6, linewidth=1.5, color='tab:blue')
ax1.axhline(d[n_latent], color='tab:red', linestyle='--', linewidth=1.5, 
            label=f'Seuil (rang effectif = {n_latent})')
ax1.fill_between(range(1, n_latent+1), d[:n_latent], alpha=0.3, color='tab:green', label='Signal')
ax1.fill_between(range(n_latent+1, len(d)+1), d[n_latent:], alpha=0.3, color='tab:orange', label='Bruit')
ax1.set_xlabel('Indice $j$')
ax1.set_ylabel('Valeur singulière $d_j$ (échelle log)')
ax1.set_title('Spectre des valeurs singulières')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim(0.5, len(d)+0.5)

# Panneau droit: variance expliquée cumulative
ax2 = axes[1]
ax2.plot(range(1, len(d)+1), cumulative_variance * 100, 'o-', markersize=6, 
         linewidth=1.5, color='tab:blue')
ax2.axhline(95, color='gray', linestyle='--', linewidth=1, label='Seuil 95%')
ax2.axvline(n_latent, color='tab:red', linestyle='--', linewidth=1.5)

# Trouver k pour 95% de variance
k_95 = np.searchsorted(cumulative_variance, 0.95) + 1
ax2.scatter([k_95], [cumulative_variance[k_95-1]*100], s=100, color='tab:red', zorder=5)
ax2.annotate(f'{k_95} composantes\npour 95%', xy=(k_95, cumulative_variance[k_95-1]*100),
             xytext=(k_95+5, cumulative_variance[k_95-1]*100-10), fontsize=10,
             arrowprops=dict(arrowstyle='->', color='gray'))

ax2.set_xlabel('Nombre de composantes $k$')
ax2.set_ylabel('Variance expliquée cumulative (%)')
ax2.set_title('Variance expliquée')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0.5, len(d)+0.5)
ax2.set_ylim(0, 105)

plt.tight_layout()
<Figure size 1200x450 with 2 Axes>

Le panneau de gauche montre le spectre des valeurs singulières en échelle logarithmique. On observe un schéma typique:

Le rang effectif est le nombre de valeurs singulières significativement au-dessus du plancher de bruit. Dans cet exemple, nous avons simulé 5 facteurs latents, et le spectre révèle bien cette structure: les 5 premières valeurs singulières dominent.

Le panneau de droite montre la variance expliquée cumulative. C’est un outil pratique pour choisir combien de composantes retenir:

De la troncature à la réduction de dimension (ACP)

L’analyse ci-dessus suggère une idée: si les dernières directions ne contiennent que du bruit, pourquoi ne pas simplement les ignorer? Au lieu de rétrécir les coefficients comme Ridge, nous pourrions tronquer la représentation en ne gardant que les kk premières directions principales.

C’est exactement l’idée de l’analyse en composantes principales (ACP). Au lieu de travailler avec les dd caractéristiques originales, nous projetons les données sur les kk premiers vecteurs singuliers:

zn=Vk(xnxˉ)Rk\mathbf{z}_n = \mathbf{V}_k^\top (\mathbf{x}_n - \bar{\mathbf{x}}) \in \mathbb{R}^k

Vk\mathbf{V}_k contient les kk premiers vecteurs singuliers (les colonnes de V\mathbf{V} correspondant aux kk plus grandes valeurs singulières).

Cette projection préserve au mieux la variance des données: les kk composantes principales capturent la direction où les données varient le plus. La reconstruction à partir de cette représentation compressée s’écrit:

x^n=Vkzn+xˉ\hat{\mathbf{x}}_n = \mathbf{V}_k \mathbf{z}_n + \bar{\mathbf{x}}

L’erreur de reconstruction est minimale parmi toutes les projections linéaires sur un sous-espace de dimension kk.

Lien entre Ridge et ACP: Les deux approches traitent le même problème (les directions à faible valeur singulière sont bruitées) mais différemment:

ApprocheTraitement des directions bruitéesType de régularisation
RidgeRétrécit (soft thresholding)Continue: garde tout, pénalise
ACPÉlimine (hard thresholding)Discrète: garde kk, ignore le reste

Ridge est appropriée pour la régression supervisée, où même les petites directions peuvent contenir du signal utile pour prédire y\mathbf{y}. L’ACP est appropriée pour la réduction de dimension non supervisée, où nous voulons une représentation compacte des données elles-mêmes.

Pourquoi λI\lambda \mathbf{I} aide

Ce terme diagonal a plusieurs effets bénéfiques:

  1. Amélioration du conditionnement: La matrice XX\mathbf{X}^\top \mathbf{X} peut être mal conditionnée (ses valeurs propres varient sur plusieurs ordres de grandeur) ou même singulière. L’ajout de λI\lambda \mathbf{I} augmente toutes les valeurs propres de λ\lambda, rendant la matrice inversible et mieux conditionnée.

  2. Rétrécissement des coefficients (shrinkage): Comme nous l’avons vu dans la section SVD ci-dessus, la solution Ridge s’écrit:

θ^ridge=j=1ddj2dj2+λujydjvj\hat{\boldsymbol{\theta}}_{\text{ridge}} = \sum_{j=1}^d \frac{d_j^2}{d_j^2 + \lambda} \frac{\mathbf{u}_j^\top \mathbf{y}}{d_j} \mathbf{v}_j

Le facteur de rétrécissement dj2dj2+λ\frac{d_j^2}{d_j^2 + \lambda} est toujours inférieur à 1, ce qui “rétrécit” chaque composante vers zéro. L’effet est différencié selon les directions:

Pour visualiser ce rétrécissement et comprendre son effet, examinons un exemple concret. L’animation suivante montre simultanément trois perspectives sur la régularisation Ridge: les données et la droite ajustée, le paysage de perte avec la contrainte, et les facteurs de rétrécissement.

Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.patches import Circle
from IPython.display import Image

# Générer des données de régression simple
np.random.seed(42)
n = 30

# Une seule caractéristique pour visualisation claire
x = np.random.uniform(-2, 2, n)
# Relation linéaire avec bruit
theta_true = 1.5
y = theta_true * x + np.random.normal(0, 0.8, n)

# Ajouter une caractéristique corrélée (pour créer de la colinéarité)
x2 = 0.9 * x + 0.3 * np.random.randn(n)

# Matrice de design avec les deux caractéristiques
X = np.column_stack([x, x2])

# Solution MCO
theta_ols = np.linalg.lstsq(X, y, rcond=None)[0]

# SVD pour analyse
U, d_svd, Vt = np.linalg.svd(X, full_matrices=False)
V = Vt.T

# Fonction pour calculer la solution Ridge
def ridge_solution(X, y, lam):
    n_features = X.shape[1]
    return np.linalg.solve(X.T @ X + lam * np.eye(n_features), X.T @ y)

# Préparer la grille pour les contours RSS
theta1_range = np.linspace(-0.5, 3, 100)
theta2_range = np.linspace(-1.5, 2, 100)
T1, T2 = np.meshgrid(theta1_range, theta2_range)

# Calculer RSS pour chaque point de la grille
RSS = np.zeros_like(T1)
for i in range(T1.shape[0]):
    for j in range(T1.shape[1]):
        theta = np.array([T1[i, j], T2[i, j]])
        residuals = y - X @ theta
        RSS[i, j] = np.sum(residuals**2)

# Créer la figure avec trois panneaux
fig = plt.figure(figsize=(15, 5))

# === Panneau 1: Données et droite ajustée ===
ax1 = fig.add_subplot(1, 3, 1)

# Données
ax1.scatter(x, y, c='tab:blue', s=50, alpha=0.7, label='Données', zorder=3)

# Grille pour tracer les droites
x_grid = np.linspace(-2.5, 2.5, 100)

# Droite MCO (fixe) - on utilise seulement theta1 car x et x2 sont très corrélés
# La prédiction effective est environ (theta1 + 0.9*theta2) * x
slope_ols = theta_ols[0] + 0.9 * theta_ols[1]  # Pente effective
y_ols = slope_ols * x_grid
ax1.plot(x_grid, y_ols, 'k-', linewidth=2, alpha=0.7, label='MCO')

# Droite Ridge (animée)
line_ridge, = ax1.plot([], [], '-', color='tab:orange', linewidth=2.5, label='Ridge')

# Ligne horizontale (prédiction = moyenne, lambda infini)
y_mean = np.mean(y)
ax1.axhline(y_mean, color='gray', linestyle=':', alpha=0.5, label=f'Moyenne ($\\lambda \\to \\infty$)')

ax1.set_xlabel('$x$')
ax1.set_ylabel('$y$')
ax1.set_title('Données et droite de régression')
ax1.legend(loc='upper left', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-2.5, 2.5)
ax1.set_ylim(-4, 5)

# Texte pour les coefficients
coef_text = ax1.text(0.98, 0.02, '', transform=ax1.transAxes, fontsize=10, 
                     ha='right', va='bottom',
                     bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

# === Panneau 2: Paysage de perte ===
ax2 = fig.add_subplot(1, 3, 2)

# Contours RSS (ellipses centrées sur OLS)
levels = np.percentile(RSS.flatten(), [5, 15, 30, 50, 70, 85, 95])
contours = ax2.contour(T1, T2, RSS, levels=levels, colors='gray', alpha=0.6)
ax2.clabel(contours, inline=True, fontsize=8, fmt='%.0f')

# Solution MCO (fixe)
ax2.plot(theta_ols[0], theta_ols[1], 'ko', markersize=12, label='MCO', zorder=5)
ax2.annotate('MCO', xy=(theta_ols[0], theta_ols[1]), 
             xytext=(theta_ols[0] + 0.2, theta_ols[1] + 0.2),
             fontsize=11, ha='left')

# Origine = coefficients nuls (prédiction constante)
ax2.plot(0, 0, 'k+', markersize=15, markeredgewidth=2, zorder=4)
ax2.annotate('$\\boldsymbol{\\theta} = 0$\n(pente nulle)', xy=(0, 0), 
             xytext=(-0.4, -1.2), fontsize=9, ha='center', color='gray')

# Cercle de contrainte Ridge (animé)
circle_ridge = Circle((0, 0), radius=np.linalg.norm(theta_ols), 
                       fill=False, edgecolor='tab:orange', linewidth=2.5, 
                       linestyle='-', alpha=0.8, zorder=3)
ax2.add_patch(circle_ridge)

# Solution Ridge (animée)
point_ridge, = ax2.plot([], [], 'o', color='tab:orange', markersize=10, 
                        label='Ridge', zorder=6)

# Chemin de régularisation
lambda_path = np.logspace(-3, 1.5, 50)
theta_path = np.array([ridge_solution(X, y, l) for l in lambda_path])
ax2.plot(theta_path[:, 0], theta_path[:, 1], 'tab:orange', linewidth=1.5, 
         alpha=0.4, linestyle='--', label='Chemin')

ax2.set_xlabel('$\\theta_1$')
ax2.set_ylabel('$\\theta_2$')
ax2.set_title('Paysage RSS et chemin de régularisation')
ax2.legend(loc='upper right', fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-0.5, 3)
ax2.set_ylim(-1.5, 2)
ax2.set_aspect('equal')

# Texte pour lambda
lambda_text = ax2.text(0.02, 0.98, '', transform=ax2.transAxes, fontsize=11, 
                       va='top', ha='left',
                       bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))

# === Panneau 3: Facteurs de rétrécissement ===
ax3 = fig.add_subplot(1, 3, 3)

lambda_range = np.linspace(0, 10, 200)
shrink1_curve = d_svd[0]**2 / (d_svd[0]**2 + lambda_range)
shrink2_curve = d_svd[1]**2 / (d_svd[1]**2 + lambda_range)

ax3.plot(lambda_range, shrink1_curve, 'b-', linewidth=2, 
         label=f'Direction forte ($d_1={d_svd[0]:.1f}$)')
ax3.plot(lambda_range, shrink2_curve, 'r-', linewidth=2, 
         label=f'Direction faible ($d_2={d_svd[1]:.2f}$)')

# Zone de surapprentissage et sous-apprentissage
ax3.axvspan(0, 0.5, alpha=0.1, color='red', label='Surapprentissage')
ax3.axvspan(5, 10, alpha=0.1, color='blue', label='Sous-apprentissage')

point_shrink1, = ax3.plot([], [], 'bo', markersize=10, zorder=3)
point_shrink2, = ax3.plot([], [], 'ro', markersize=10, zorder=3)

ax3.axhline(1.0, color='gray', linestyle='--', alpha=0.5)
ax3.set_xlabel('$\\lambda$')
ax3.set_ylabel('Facteur de rétrécissement')
ax3.set_title('Rétrécissement par direction SVD')
ax3.legend(loc='center right', fontsize=8)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 10)
ax3.set_ylim(0, 1.1)

shrink_text = ax3.text(0.02, 0.5, '', transform=ax3.transAxes, fontsize=10, 
                       va='center', ha='left',
                       bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()

# Fonction d'animation
def animate(frame):
    if frame < 80:
        lam = (frame / 80) * 10
    else:
        lam = 10
    
    # Solution Ridge
    theta_ridge = ridge_solution(X, y, lam)
    
    # Panneau 1: Mettre à jour la droite
    slope_ridge = theta_ridge[0] + 0.9 * theta_ridge[1]
    y_ridge = slope_ridge * x_grid
    line_ridge.set_data(x_grid, y_ridge)
    coef_text.set_text(f'Pente MCO: {slope_ols:.2f}\nPente Ridge: {slope_ridge:.2f}')
    
    # Panneau 2: Mettre à jour le cercle et le point
    norm_ridge = np.linalg.norm(theta_ridge)
    circle_ridge.set_radius(norm_ridge)
    point_ridge.set_data([theta_ridge[0]], [theta_ridge[1]])
    lambda_text.set_text(f'$\\lambda = {lam:.1f}$')
    
    # Panneau 3: Mettre à jour les points de rétrécissement
    shrink1 = d_svd[0]**2 / (d_svd[0]**2 + lam)
    shrink2 = d_svd[1]**2 / (d_svd[1]**2 + lam)
    point_shrink1.set_data([lam], [shrink1])
    point_shrink2.set_data([lam], [shrink2])
    shrink_text.set_text(f'Facteur dir. 1: {shrink1:.2f}\nFacteur dir. 2: {shrink2:.2f}')
    
    return (line_ridge, point_ridge, circle_ridge, lambda_text, 
            point_shrink1, point_shrink2, coef_text, shrink_text)

# Créer l'animation
anim = FuncAnimation(fig, animate, frames=90, interval=80, blit=False, repeat=True)
anim.save('_static/ridge_geometry.gif', writer='pillow', fps=12, dpi=100)
plt.close()

# Afficher le GIF
Image(filename='_static/ridge_geometry.gif')
<IPython.core.display.Image object>

L’animation relie trois perspectives sur la régularisation Ridge lorsque λ\lambda augmente de 0 à 10:

Panneau de gauche (données et ajustement): Les points bleus sont les données d’entraînement. La droite noire est l’ajustement MCO (λ=0\lambda = 0), la droite orange est l’ajustement Ridge. À mesure que λ\lambda augmente, la pente de la droite Ridge diminue, se rapprochant de la ligne horizontale (prédiction constante égale à la moyenne). C’est le rétrécissement vers zéro: Ridge “tire” les coefficients vers l’origine, ce qui réduit la pente.

Panneau central (paysage de perte): Chaque point de ce plan représente un choix de coefficients (θ1,θ2)(\theta_1, \theta_2). Les contours gris montrent la fonction de coût RSS: plus on est proche du point noir (MCO), plus l’erreur sur les données d’entraînement est faible. L’ellipse est allongée car x1x_1 et x2x_2 sont corrélées (colinéarité). L’origine θ=(0,0)\boldsymbol{\theta} = (0, 0) correspond à une pente nulle (prédiction constante). Le cercle orange montre la norme de la solution Ridge courante θ^ridge2\|\hat{\boldsymbol{\theta}}_{\text{ridge}}\|_2. La formulation pénalisée RSS+λθ2\text{RSS} + \lambda\|\boldsymbol{\theta}\|^2 est équivalente à la formulation contrainte minRSS\min \text{RSS} sous θ2t\|\boldsymbol{\theta}\|^2 \leq t, où λ\lambda joue le rôle du multiplicateur de Lagrange: pour chaque λ\lambda, il existe un tt tel que les deux problèmes ont la même solution. À mesure que λ\lambda augmente, la solution se déplace le long du chemin de régularisation vers l’origine.

Panneau de droite (rétrécissement différencié): La direction “forte” (grande valeur singulière d1d_1, où les données sont dispersées) est peu affectée par la régularisation. La direction “faible” (petite valeur singulière d2d_2, direction de colinéarité) est rétrécit beaucoup plus rapidement. C’est le cœur de l’effet Ridge: pénaliser davantage les directions où le signal est faible et l’estimation instable.

L’intuition géométrique est la suivante: quand les données sont colinéaires, l’ellipse RSS est très allongée. De petites perturbations dans les données causent de grands déplacements de la solution MCO le long de l’axe allongé. La pénalité Ridge ajoute un terme λθ2\lambda\|\boldsymbol{\theta}\|^2 qui « tire » la solution vers l’origine. Dans la formulation contrainte équivalente, cela correspond à chercher le minimum de RSS à l’intérieur d’une boule de rayon t\sqrt{t}. Plus la boule est petite (plus λ\lambda est grand), plus la solution est proche de l’origine et donc plus stable.

  1. Stabilité numérique: Quand XX\mathbf{X}^\top \mathbf{X} est presque singulière, de petites perturbations dans les données causent de grandes variations dans θ^MCO\hat{\boldsymbol{\theta}}_{\text{MCO}}. La régularisation réduit cette sensibilité.

Classes de modèles et expansion de caractéristiques

Les exemples de régularisation ci-dessus utilisaient des caractéristiques polynomiales: au lieu de prédire yy directement à partir de xx, nous avons construit des caractéristiques [1,x,x2,,x15][1, x, x^2, \ldots, x^{15}] et appliqué un modèle linéaire dans cet espace étendu. Cette technique s’appelle l’expansion de caractéristiques et mérite d’être formalisée.

Trois familles de modèles

Situons les modèles linéaires dans une hiérarchie plus large. Nous distinguons trois familles de complexité croissante:

  1. Modèles linéaires: f(x;θ)=θx+bf(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \mathbf{x} + b. La sortie est une combinaison linéaire des entrées. Simple, interprétable, mais limité aux relations linéaires.

  2. Modèles à expansion de caractéristiques: f(x;θ)=θϕ(x)+bf(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \boldsymbol{\phi}(\mathbf{x}) + b, où ϕ:RdRD\boldsymbol{\phi}: \mathbb{R}^d \to \mathbb{R}^D est une transformation non linéaire fixée à l’avance (par exemple, polynomiale). Le modèle reste linéaire dans les paramètres θ\boldsymbol{\theta}, ce qui facilite l’optimisation, mais peut capturer des relations non linéaires en x\mathbf{x}. L’espace de redescription a souvent une dimension DdD \gg d.

  3. Réseaux de neurones: f(x;θ)=fK(fK1(f1(x;θ1);θK1);θK)f(\mathbf{x}; \boldsymbol{\theta}) = f_K(f_{K-1}(\cdots f_1(\mathbf{x}; \boldsymbol{\theta}_1); \boldsymbol{\theta}_{K-1}); \boldsymbol{\theta}_K). Une composition de KK fonctions non linéaires, chacune avec ses propres paramètres. Contrairement aux modèles à expansion fixe, les réseaux de neurones apprennent la représentation ϕ\boldsymbol{\phi} en même temps que les paramètres θ\boldsymbol{\theta}.

Cette progression capture l’évolution historique du domaine: des modèles linéaires classiques aux méthodes à noyaux (expansion implicite), puis aux réseaux profonds qui apprennent leurs propres représentations. Nous verrons les réseaux de neurones en détail dans les chapitres suivants; concentrons-nous ici sur les deux premières familles.

Expansion de caractéristiques

Pour capturer des relations non linéaires tout en gardant un modèle linéaire dans les paramètres, nous transformons les entrées. En régression polynomiale, nous appliquons une fonction ϕ:RRk+1\phi: \mathbb{R} \to \mathbb{R}^{k+1}:

ϕ(x)=[1,x,x2,,xk]\phi(x) = [1, x, x^2, \ldots, x^k]

La prédiction devient f(x;θ)=θϕ(x)f(x; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \phi(x). Le modèle est polynomial en xx mais linéaire en θ\boldsymbol{\theta}, ce qui permet d’utiliser les mêmes algorithmes d’optimisation (MCO, Ridge).

Le degré kk contrôle la capacité du modèle. Avec k=1k = 1, nous avons une droite. Avec kk élevé, le polynôme peut osciller pour passer par tous les points d’entraînement. Avec k=N1k = N - 1, nous pouvons interpoler exactement les NN points: le risque empirique atteint zéro. Mais un polynôme qui passe exactement par les points d’entraînement n’a aucune raison de bien prédire les nouveaux points.

Illustrons ce phénomène avec les données de freinage. Nous ajustons des polynômes de degrés 1, 2, 5 et 15, et comparons leurs erreurs sur les ensembles d’entraînement et de test.

Source
import numpy as np
import matplotlib.pyplot as plt
import warnings

# Suppress polyfit warnings for high-degree polynomials (expected for this demo)
warnings.filterwarnings('ignore', message='Polyfit may be poorly conditioned')

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

# Train/test split
np.random.seed(42)
indices = np.random.permutation(len(speed))
train_idx, test_idx = indices[:35], indices[35:]
speed_train, dist_train = speed[train_idx], dist[train_idx]
speed_test, dist_test = speed[test_idx], dist[test_idx]

degrees_to_plot = [1, 2, 5, 15]
degrees_eval = range(1, 16)
fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# Pre-compute all errors for the summary plot later
all_train_errors = []
all_test_errors = []
for deg in degrees_eval:
    coeffs = np.polyfit(speed_train, dist_train, deg)
    all_train_errors.append(np.mean((dist_train - np.polyval(coeffs, speed_train))**2))
    all_test_errors.append(np.mean((dist_test - np.polyval(coeffs, speed_test))**2))

for ax, deg in zip(axes.flat, degrees_to_plot):
    # Fit polynomial
    coeffs = np.polyfit(speed_train, dist_train, deg)
    
    # Predictions
    pred_train = np.polyval(coeffs, speed_train)
    pred_test = np.polyval(coeffs, speed_test)
    
    # MSE
    mse_train = np.mean((dist_train - pred_train)**2)
    mse_test = np.mean((dist_test - pred_test)**2)
    
    # Plot
    ax.scatter(speed_train, dist_train, alpha=0.6, s=30, label='Entraînement')
    ax.scatter(speed_test, dist_test, alpha=0.6, s=30, marker='s', label='Test')
    
    speed_grid = np.linspace(3, 26, 200)
    pred_grid = np.polyval(coeffs, speed_grid)
    # Clip extreme predictions for visualization
    pred_grid = np.clip(pred_grid, -50, 200)
    ax.plot(speed_grid, pred_grid, 'k-', alpha=0.7)
    
    ax.set_xlim(3, 26)
    ax.set_ylim(-20, 150)
    ax.set_xlabel('Vitesse (mph)')
    ax.set_ylabel('Distance (ft)')
    ax.set_title(f'Degré {deg}: Entr. EQM={mse_train:.1f}, Test EQM={mse_test:.1f}')
    if deg == 1:
        ax.legend()

plt.tight_layout()
<Figure size 1000x800 with 4 Axes>

Le polynôme de degré 1 (droite) ne capture pas la courbure des données: c’est du sous-apprentissage. Le polynôme de degré 2 capture bien la relation quadratique. Le polynôme de degré 5 commence à osciller. Le polynôme de degré 15 passe près de tous les points d’entraînement, mais ses oscillations produisent des prédictions absurdes entre les points: c’est du surapprentissage.

Source
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(degrees_eval, all_train_errors, 'o-', linewidth=2, label='Erreur entraînement')
ax.plot(degrees_eval, all_test_errors, 's-', linewidth=2, label='Erreur test')

# Utiliser une échelle logarithmique car l'erreur de test explose
ax.set_yscale('log')

ax.set_xlabel('Degré du polynôme (complexité)')
ax.set_ylabel('EQM (échelle log)')
ax.set_xticks(range(1, 16, 2))
ax.grid(True, which="both", ls="-", alpha=0.2)
ax.legend()

ax.set_title('Compromis biais-variance')
plt.tight_layout()
<Figure size 800x500 with 1 Axes>

L’erreur d’entraînement diminue avec le degré du polynôme. L’erreur de test diminue d’abord (quand le modèle gagne en expressivité), puis augmente (quand le modèle commence à mémoriser le bruit). Le meilleur modèle se trouve à l’intersection de ces deux tendances. La régularisation Ridge, vue précédemment, est une alternative au choix du degré: elle permet d’utiliser un modèle de haute capacité tout en contrôlant le surapprentissage.

Décomposition biais-variance

Ce compromis peut être formalisé mathématiquement. Supposons que les données suivent le modèle y=f(x)+ϵy = f^*(\mathbf{x}) + \epsilon, où ff^* est la vraie fonction (le prédicteur de Bayes optimal pour la perte quadratique), et ϵ\epsilon est un bruit de moyenne nulle et de variance σ2\sigma^2.

Clarification sur les variables aléatoires: Dans cette analyse, il y a trois sources d’aléa qu’il faut bien distinguer:

Autrement dit, l’espérance E[f^(x)]\mathbb{E}[\hat{f}(\mathbf{x})] moyenne sur tous les échantillons d’entraînement possibles: si nous pouvions répéter l’expérience d’apprentissage un grand nombre de fois avec différents D\mathcal{D}, quelle serait la prédiction moyenne pour ce x\mathbf{x}?

L’erreur quadratique moyenne, moyennée sur les échantillons d’entraînement possibles et le bruit, se décompose comme suit:

ED,ϵ[(f^(x)y)2]=E[(f^(x)f(x)ϵ)2]\mathbb{E}_{\mathcal{D}, \epsilon}[(\hat{f}(\mathbf{x}) - y)^2] = \mathbb{E}[(\hat{f}(\mathbf{x}) - f^*(\mathbf{x}) - \epsilon)^2]

En développant le carré et utilisant E[ϵ]=0\mathbb{E}[\epsilon] = 0 ainsi que l’indépendance entre ϵ\epsilon et f^\hat{f}:

=E[(f^(x)f(x))2]+σ2= \mathbb{E}[(\hat{f}(\mathbf{x}) - f^*(\mathbf{x}))^2] + \sigma^2

Pour le premier terme, ajoutons et retranchons E[f^(x)]\mathbb{E}[\hat{f}(\mathbf{x})]:

E[(f^(x)f(x))2]=E[(f^(x)E[f^(x)]+E[f^(x)]f(x))2]\mathbb{E}[(\hat{f}(\mathbf{x}) - f^*(\mathbf{x}))^2] = \mathbb{E}[(\hat{f}(\mathbf{x}) - \mathbb{E}[\hat{f}(\mathbf{x})] + \mathbb{E}[\hat{f}(\mathbf{x})] - f^*(\mathbf{x}))^2]

En développant et utilisant E[f^(x)E[f^(x)]]=0\mathbb{E}[\hat{f}(\mathbf{x}) - \mathbb{E}[\hat{f}(\mathbf{x})]] = 0:

=E[(f^(x)E[f^(x)])2]Var(f^(x))+(E[f^(x)]f(x))2Biais2(f^(x))= \underbrace{\mathbb{E}[(\hat{f}(\mathbf{x}) - \mathbb{E}[\hat{f}(\mathbf{x})])^2]}_{\text{Var}(\hat{f}(\mathbf{x}))} + \underbrace{(\mathbb{E}[\hat{f}(\mathbf{x})] - f^*(\mathbf{x}))^2}_{\text{Biais}^2(\hat{f}(\mathbf{x}))}

Nous obtenons la décomposition biais-variance:

E[(f^(x)y)2]=Biais2(f^(x))+Var(f^(x))+σ2\boxed{\mathbb{E}[(\hat{f}(\mathbf{x}) - y)^2] = \text{Biais}^2(\hat{f}(\mathbf{x})) + \text{Var}(\hat{f}(\mathbf{x})) + \sigma^2}

Chaque terme a une interprétation précise:

Cette décomposition explique pourquoi l’erreur de test a une forme en U: à faible complexité, le biais domine; à haute complexité, la variance domine. Le minimum se trouve au point où la somme des deux est minimale.

Intuition géométrique: pourquoi la dimension supérieure aide

L’expansion de caractéristiques semble être un simple changement de variables, mais elle cache une idée géométrique profonde. Pour comprendre pourquoi projeter les données dans un espace de dimension supérieure permet de capturer des relations non linéaires, examinons d’abord le cas de la régression, puis celui de la classification.

Le plan caché derrière la parabole

Considérons une régression quadratique: f(x)=θ0+θ1x+θ2x2f(x) = \theta_0 + \theta_1 x + \theta_2 x^2. Cette fonction est non linéaire en xx (c’est une parabole), mais linéaire dans les paramètres θ=(θ0,θ1,θ2)\boldsymbol{\theta} = (\theta_0, \theta_1, \theta_2). Que signifie cette distinction géométriquement?

Introduisons l’espace des caractéristiques ϕ(x)=(1,x,x2)\phi(x) = (1, x, x^2). Chaque valeur de xx correspond à un point dans R3\mathbb{R}^3. Ces points ne sont pas dispersés arbitrairement: ils vivent sur une courbe particulière, la courbe des moments (moment curve), qui ressemble à une rampe tordue.

Source
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Créer la courbe des moments: (x, x², x³) pour visualisation
# On utilise (1, x, x²) mais on projette sur (x, x², y) pour la visualisation
fig = plt.figure(figsize=(12, 5))

# Gauche: les fonctions de base
ax1 = fig.add_subplot(121)
x = np.linspace(-2, 2, 100)
ax1.plot(x, np.ones_like(x), 'b-', linewidth=2, label=r'$\phi_0(x) = 1$')
ax1.plot(x, x, 'orange', linewidth=2, label=r'$\phi_1(x) = x$')
ax1.plot(x, x**2, 'g-', linewidth=2, label=r'$\phi_2(x) = x^2$')
ax1.axhline(0, color='gray', linewidth=0.5)
ax1.axvline(0, color='gray', linewidth=0.5)
ax1.set_xlabel('$x$')
ax1.set_ylabel(r'$\phi_j(x)$')
ax1.set_title('Les fonctions de base')
ax1.legend()
ax1.set_ylim(-2.5, 4.5)
ax1.grid(True, alpha=0.3)

# Droite: combinaisons linéaires
ax2 = fig.add_subplot(122)
x = np.linspace(-2, 2, 100)

# Différentes combinaisons de coefficients
combinations = [
    ((1, 0, 0), 'Constante: $1$'),
    ((0, 1, 0), 'Linéaire: $x$'),
    ((0, 0, 1), 'Quadratique: $x^2$'),
    ((1, -0.5, 0.5), 'Combinaison: $1 - 0.5x + 0.5x^2$'),
]

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red']
for (theta0, theta1, theta2), label in combinations:
    y = theta0 + theta1 * x + theta2 * x**2
    ax2.plot(x, y, linewidth=2, label=label)

ax2.axhline(0, color='gray', linewidth=0.5)
ax2.axvline(0, color='gray', linewidth=0.5)
ax2.set_xlabel('$x$')
ax2.set_ylabel('$f(x)$')
ax2.set_title('Combinaisons linéaires des fonctions de base')
ax2.legend(loc='upper center')
ax2.set_ylim(-2.5, 4.5)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
<Figure size 1200x500 with 2 Axes>

Les fonctions de base {1,x,x2}\{1, x, x^2\} sont les “ingrédients” du modèle. La régression polynomiale cherche les coefficients θ0,θ1,θ2\theta_0, \theta_1, \theta_2 qui mélangent ces ingrédients de façon optimale. Chaque combinaison produit une courbe différente, mais toutes sont des paraboles (ou des cas dégénérés: droites, constantes).

Voici l’insight géométrique clé: dans l’espace (x,x2,y)(x, x^2, y), le modèle y=θ0+θ1x+θ2x2y = \theta_0 + \theta_1 x + \theta_2 x^2 définit un plan. La parabole que nous voyons dans le graphique (x,y)(x, y) est simplement la projection de ce plan sur notre espace de visualisation.

Source
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import Image

# Générer des données avec relation quadratique
np.random.seed(42)
n_points = 40
x_data = np.random.uniform(-1.8, 1.8, n_points)
y_true = 0.5 + 0.3 * x_data + 0.8 * x_data**2
y_data = y_true + np.random.normal(0, 0.3, n_points)

# Ajuster le modèle quadratique
X_design = np.column_stack([np.ones(n_points), x_data, x_data**2])
theta = np.linalg.lstsq(X_design, y_data, rcond=None)[0]

# Créer l'animation
fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111, projection='3d')

# Grille pour le plan de régression
x_grid = np.linspace(-2, 2, 20)
x2_grid = np.linspace(0, 4, 20)
X_plane, X2_plane = np.meshgrid(x_grid, x2_grid)
Y_plane = theta[0] + theta[1] * X_plane + theta[2] * X2_plane

# Surface parabolique z = x²
x_surf = np.linspace(-2, 2, 30)
X_surf, Y_surf_temp = np.meshgrid(x_surf, np.linspace(-1, 5, 30))
Z_surf = X_surf**2

def init():
    ax.clear()
    return []

def animate(frame):
    ax.clear()
    
    # Animation en trois phases
    if frame < 15:
        # Phase 1: Vue 2D (de côté, cachant x²)
        elev = 0
        azim = 0
        show_plane = False
        show_surface = False
        title = 'Vue 2D: régression quadratique'
    elif frame < 35:
        # Phase 2: Rotation révélant la 3ème dimension
        progress = (frame - 15) / 20
        elev = progress * 25
        azim = progress * 45
        show_plane = False
        show_surface = True
        title = 'Rotation: découverte de la dimension $x^2$...'
    else:
        # Phase 3: Vue 3D avec le plan
        elev = 25
        azim = 45 + (frame - 35) * 2
        show_plane = True
        show_surface = True
        title = r'Le plan $y = \theta_0 + \theta_1 x + \theta_2 x^2$ dans lespace 3D'
    
    ax.view_init(elev=elev, azim=azim)
    
    # Surface parabolique (rampe z = x²)
    if show_surface:
        ax.plot_surface(X_surf, Z_surf, Y_surf_temp, alpha=0.1, color='gray')
    
    # Points de données dans l'espace 3D: (x, x², y)
    ax.scatter(x_data, x_data**2, y_data, c='tab:blue', s=50, alpha=0.8, 
               label='Données', depthshade=True)
    
    # Plan de régression
    if show_plane:
        ax.plot_surface(X_plane, X2_plane, Y_plane, alpha=0.4, color='tab:orange',
                       label='Plan de régression')
    
    # Courbe de régression sur la surface z = x²
    x_curve = np.linspace(-1.8, 1.8, 100)
    y_curve = theta[0] + theta[1] * x_curve + theta[2] * x_curve**2
    ax.plot(x_curve, x_curve**2, y_curve, 'r-', linewidth=3, 
            label='Courbe ajustée')
    
    ax.set_xlabel('$x$')
    ax.set_ylabel('$x^2$')
    ax.set_zlabel('$y$')
    ax.set_xlim(-2, 2)
    ax.set_ylim(0, 4)
    ax.set_zlim(-1, 5)
    ax.set_title(title, fontsize=11)
    
    return []

anim = FuncAnimation(fig, animate, init_func=init, frames=70, interval=100, blit=True)
anim.save('_static/regression_plane_3d.gif', writer='pillow', fps=10, dpi=100)
plt.close()

# Afficher le GIF
Image(filename='_static/regression_plane_3d.gif')
<IPython.core.display.Image object>

L’animation révèle la structure cachée de la régression polynomiale:

  1. Vue 2D initiale: On voit les données et la parabole ajustée, comme dans un graphique classique.

  2. Rotation: En faisant pivoter la vue, on découvre que chaque point (xi,yi)(x_i, y_i) vit en réalité dans un espace 3D aux coordonnées (xi,xi2,yi)(x_i, x_i^2, y_i). La surface grise représente la “rampe” z=x2z = x^2 sur laquelle tous les points sont contraints de vivre.

  3. Le plan de régression: Le modèle y=θ0+θ1x+θ2x2y = \theta_0 + \theta_1 x + \theta_2 x^2 est un plan (en orange) dans cet espace 3D. Ce plan est choisi pour minimiser les distances verticales aux points.

  4. La courbe ajustée: La parabole rouge est l’intersection du plan avec la rampe z=x2z = x^2. C’est ce que nous voyons quand nous projetons le tout sur le plan (x,y)(x, y).

Cette perspective unifie le “linéaire dans les paramètres” et le “non linéaire en xx”:

De la régression à la classification

La même intuition s’applique en classification, avec une différence: au lieu de chercher un plan qui ajuste les données, on cherche un hyperplan qui les sépare. Voyons comment l’expansion de caractéristiques transforme des données non linéairement séparables en données linéairement séparables.

De 1D à 2D: séparer l’inséparable

Considérons des points sur une droite, répartis en deux classes: les points bleus au centre, les points orange aux extrémités. Aucun seuil unique ne peut séparer ces deux classes.

Source
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# Classe bleue: points au centre
x_blue = np.random.uniform(-0.5, 0.5, 15)
# Classe orange: points aux extrémités
x_orange = np.concatenate([np.random.uniform(-1.5, -0.8, 8), 
                           np.random.uniform(0.8, 1.5, 8)])

fig, ax = plt.subplots(figsize=(10, 2))
ax.scatter(x_blue, np.zeros_like(x_blue), c='tab:blue', s=80, label='Classe A', zorder=3)
ax.scatter(x_orange, np.zeros_like(x_orange), c='tab:orange', s=80, label='Classe B', zorder=3)
ax.axhline(0, color='gray', linewidth=0.5, zorder=1)
ax.set_xlim(-2, 2)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('$x$')
ax.set_yticks([])
ax.legend(loc='upper right')
ax.set_title('Données 1D: aucun seuil ne sépare les deux classes')
plt.tight_layout()
<Figure size 1000x200 with 1 Axes>

Appliquons maintenant l’expansion ϕ(x)=(x,x2)\phi(x) = (x, x^2). Chaque point est projeté sur une parabole dans l’espace 2D. Les points proches de zéro (classe bleue) restent bas, tandis que les points éloignés de zéro (classe orange) montent.

Source
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Gauche: espace original avec tentative de séparation
ax = axes[0]
ax.scatter(x_blue, np.zeros_like(x_blue), c='tab:blue', s=80, zorder=3)
ax.scatter(x_orange, np.zeros_like(x_orange), c='tab:orange', s=80, zorder=3)
ax.axhline(0, color='gray', linewidth=0.5, zorder=1)
ax.axvline(0.3, color='red', linestyle='--', linewidth=2, label='Seuil?')
ax.set_xlim(-2, 2)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('$x$')
ax.set_yticks([])
ax.set_title('Espace original: pas de séparation linéaire')
ax.legend()

# Droite: espace transformé
ax = axes[1]
ax.scatter(x_blue, x_blue**2, c='tab:blue', s=80, zorder=3, label='Classe A')
ax.scatter(x_orange, x_orange**2, c='tab:orange', s=80, zorder=3, label='Classe B')

# Ligne de séparation dans l'espace transformé
x_line = np.linspace(-2, 2, 100)
threshold = 0.6
ax.axhline(threshold, color='green', linestyle='-', linewidth=2, label='Frontière linéaire')
ax.fill_between(x_line, 0, threshold, alpha=0.1, color='blue')
ax.fill_between(x_line, threshold, 2.5, alpha=0.1, color='orange')

# Parabole de référence
x_curve = np.linspace(-1.6, 1.6, 100)
ax.plot(x_curve, x_curve**2, 'k--', alpha=0.3, linewidth=1)

ax.set_xlim(-2, 2)
ax.set_ylim(-0.1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$x^2$')
ax.set_title(r'Espace transformé $\phi(x) = (x, x^2)$: séparation linéaire!')
ax.legend()

plt.tight_layout()
<Figure size 1200x400 with 2 Axes>

Dans l’espace transformé, une simple droite horizontale sépare les deux classes. Cette droite correspond, dans l’espace original, à deux seuils: x2<0.6x^2 < 0.6, soit x<0.60.77|x| < \sqrt{0.6} \approx 0.77. L’expansion de caractéristiques a transformé une frontière de décision non linéaire (un intervalle) en une frontière linéaire (une droite).

De 2D à 3D: soulever pour séparer

Passons à un exemple plus visuel. Considérons deux classes disposées en cercles concentriques: la classe bleue forme un disque central, la classe orange forme un anneau extérieur. Aucune droite ne peut séparer ces deux régions.

Source
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Classe bleue: disque central
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)

# Classe orange: anneau extérieur
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x_blue, y_blue, c='tab:blue', s=40, alpha=0.7, label='Classe A')
ax.scatter(x_orange, y_orange, c='tab:orange', s=40, alpha=0.7, label='Classe B')

# Montrer qu'une droite ne peut pas séparer
theta_line = np.pi/4
x_line = np.linspace(-2, 2, 100)
y_line = np.tan(theta_line) * x_line
ax.plot(x_line, y_line, 'r--', linewidth=2, alpha=0.7, label='Droite?')

ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_aspect('equal')
ax.legend()
ax.set_title('Cercles concentriques: pas de séparation linéaire en 2D')
plt.tight_layout()
<Figure size 600x600 with 1 Axes>

Appliquons l’expansion ϕ(x1,x2)=(x1,x2,x12+x22)\phi(x_1, x_2) = (x_1, x_2, x_1^2 + x_2^2). La troisième coordonnée est le carré de la distance à l’origine: z=r2z = r^2. Les points proches du centre (petit rr) sont “soulevés” moins haut que les points éloignés (grand rr).

L’animation suivante montre cette transformation. En faisant pivoter la vue, on voit que les données, une fois projetées dans l’espace 3D, deviennent séparables par un plan horizontal.

Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import Image

np.random.seed(123)

# Classe bleue: disque central
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)
z_blue = x_blue**2 + y_blue**2

# Classe orange: anneau extérieur  
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)
z_orange = x_orange**2 + y_orange**2

# Créer la figure
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

# Plan de séparation
z_sep = 0.75
xx, yy = np.meshgrid(np.linspace(-1.8, 1.8, 10), np.linspace(-1.8, 1.8, 10))
zz = np.ones_like(xx) * z_sep

def init():
    ax.clear()
    return []

def animate(frame):
    ax.clear()
    
    # Élévation: commence à 90° (vue de dessus), descend à 25°
    if frame < 20:
        elev = 90 - frame * 3.25  # 90 -> 25
        azim = 45
    else:
        elev = 25
        azim = 45 + (frame - 20) * 4  # rotation horizontale
    
    ax.view_init(elev=elev, azim=azim)
    
    # Points
    ax.scatter(x_blue, y_blue, z_blue, c='tab:blue', s=30, alpha=0.8, label='Classe A')
    ax.scatter(x_orange, y_orange, z_orange, c='tab:orange', s=30, alpha=0.8, label='Classe B')
    
    # Plan de séparation (apparaît après la descente)
    if frame >= 15:
        alpha_plane = min(0.3, (frame - 15) * 0.06)
        ax.plot_surface(xx, yy, zz, alpha=alpha_plane, color='green')
    
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$x_1^2 + x_2^2$')
    ax.set_xlim(-1.8, 1.8)
    ax.set_ylim(-1.8, 1.8)
    ax.set_zlim(0, 2.5)
    
    # Titre dynamique
    if frame < 10:
        ax.set_title('Vue de dessus: cercles concentriques', fontsize=11)
    elif frame < 20:
        ax.set_title('Rotation: on découvre la structure 3D...', fontsize=11)
    else:
        ax.set_title(r'Espace $\phi(x_1,x_2) = (x_1, x_2, x_1^2+x_2^2)$: un plan sépare!', fontsize=11)
    
    return []

anim = FuncAnimation(fig, animate, init_func=init, frames=60, interval=100, blit=True)
anim.save('_static/feature_expansion_3d.gif', writer='pillow', fps=10, dpi=100)
plt.close()

# Afficher le GIF
Image(filename='_static/feature_expansion_3d.gif')
<IPython.core.display.Image object>

L’animation montre comment l’expansion de caractéristiques “déplie” la géométrie des données. Vue de dessus, la structure est celle des cercles concentriques originaux. Mais la troisième dimension z=x12+x22z = x_1^2 + x_2^2 sépare verticalement les deux classes: le disque central reste bas, l’anneau extérieur monte. Un plan horizontal (en vert) suffit alors à séparer les classes.

Ce plan horizontal z=0.75z = 0.75 correspond, dans l’espace original 2D, à un cercle de rayon 0.750.87\sqrt{0.75} \approx 0.87. La frontière de décision linéaire en 3D devient une frontière circulaire en 2D.

Source
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Recréer les données (même seed que l'animation)
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)

n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x_blue, y_blue, c='tab:blue', s=40, alpha=0.7, label='Classe A')
ax.scatter(x_orange, y_orange, c='tab:orange', s=40, alpha=0.7, label='Classe B')

# Cercle de décision (projection du plan z = 0.75)
theta_circle = np.linspace(0, 2*np.pi, 100)
r_decision = np.sqrt(0.75)
ax.plot(r_decision * np.cos(theta_circle), r_decision * np.sin(theta_circle), 
        'g-', linewidth=2.5, label=f'Frontière: $r = {r_decision:.2f}$')

ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_aspect('equal')
ax.legend()
ax.set_title('Frontière de décision projetée en 2D')
plt.tight_layout()
<Figure size 600x600 with 1 Axes>

Le principe unificateur

Ces exemples de régression et de classification illustrent le même principe géométrique:

RégressionClassification
ObjectifAjuster les donnéesSéparer les classes
Dans l’espace originalCourbe (parabole, etc.)Frontière courbe (cercle, etc.)
Dans l’espace des caractéristiquesHyperplan d’ajustementHyperplan séparateur
La courbe/frontière est...L’intersection du plan avec la surface ϕ(x)\phi(x)La projection de l’hyperplan

L’expansion de caractéristiques transforme un problème non linéaire en un problème linéaire dans un espace de dimension supérieure. Les modèles linéaires, simples à optimiser et à analyser, deviennent alors suffisants pour capturer des structures complexes.

En augmentant la dimension de l’espace de représentation, nous augmentons la capacité du modèle. Pour la classification, un résultat classique de géométrie affirme que NN points en position générale sont presque sûrement séparables par un hyperplan si la dimension de l’espace est au moins NN. Pour la régression, un polynôme de degré N1N-1 peut interpoler exactement NN points.

Mais cette flexibilité a un coût: plus l’espace est grand, plus le modèle risque de mémoriser les particularités des données d’entraînement plutôt que d’apprendre la structure sous-jacente. C’est le compromis biais-variance, et c’est pourquoi la régularisation (vue précédemment) est si importante pour les modèles à haute capacité.

Évaluation et choix de modèle

En pratique, nous estimons le risque par le risque empirique sur un ensemble de test Dtest\mathcal{D}_{\text{test}} disjoint de l’ensemble d’entraînement. Un troisième ensemble, l’ensemble de validation, sert à choisir parmi plusieurs modèles ou à régler des hyperparamètres. L’ensemble de test doit rester intact jusqu’à l’évaluation finale, pour fournir une estimation non biaisée.

Cette séparation est importante. Si nous utilisons l’ensemble de test pour faire des choix (quel modèle garder, quelle valeur d’hyperparamètre utiliser), l’estimation de performance sur ce même ensemble devient optimiste.

Hyperparamètres et validation

De nombreux modèles ont des hyperparamètres: des choix qui doivent être faits avant l’entraînement et qui ne sont pas appris à partir des données. Le degré kk d’un polynôme, le nombre de voisins dans les kk plus proches voisins, ou le coefficient de régularisation λ\lambda sont des exemples d’hyperparamètres.

Un hyperparamètre mal choisi peut mener au surapprentissage (modèle trop complexe) ou au sous-apprentissage (modèle trop simple). La méthode standard pour choisir un hyperparamètre est la validation: on réserve une partie des données (typiquement 20%) comme ensemble de validation, on entraîne le modèle pour plusieurs valeurs de l’hyperparamètre, et on retient celle qui minimise l’erreur sur l’ensemble de validation.

Plus formellement, pour un hyperparamètre hh, définissons le risque de validation:

R^hval=R^(f^h(Dtrain),Dvalid)\hat{\mathcal{R}}^{\text{val}}_h = \hat{\mathcal{R}}\left(\hat{f}_h(\mathcal{D}_{\text{train}}), \mathcal{D}_{\text{valid}}\right)

f^h(Dtrain)\hat{f}_h(\mathcal{D}_{\text{train}}) est le modèle entraîné avec l’hyperparamètre hh. La recherche par grille consiste à évaluer ce risque pour un ensemble de valeurs candidates et à retenir:

h=argminh{h1,,hK}R^hvalh^* = \arg\min_{h \in \{h_1, \ldots, h_K\}} \hat{\mathcal{R}}^{\text{val}}_h

Une fois hh^* choisi, on peut ré-entraîner le modèle sur l’ensemble des données (entraînement + validation) pour obtenir le modèle final.

Validation croisée

Quand les données sont peu nombreuses, réserver 20% pour la validation peut être coûteux. La validation croisée offre une alternative.

L’idée est de partitionner les données en KK blocs. Pour chaque bloc kk, on entraîne le modèle sur les K1K-1 autres blocs et on évalue sur le bloc kk. Le risque de validation croisée est la moyenne des KK évaluations:

R^hcv=1Kk=1KR^(f^h(Dk),Dk)\hat{\mathcal{R}}^{\text{cv}}_h = \frac{1}{K} \sum_{k=1}^K \hat{\mathcal{R}}\left(\hat{f}_h(\mathcal{D}_{-k}), \mathcal{D}_k\right)

Dk\mathcal{D}_{-k} désigne toutes les données sauf le bloc kk.

Source
import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 3))

K = 5
for fold in range(K):
    for k in range(K):
        if k == fold:
            ax.barh(fold, 1, left=k, color='C1', edgecolor='black', linewidth=1, label='Validation' if fold == 0 else '')
        else:
            ax.barh(fold, 1, left=k, color='C0', edgecolor='black', linewidth=1, label='Entraînement' if fold == 0 and k == 0 else '')

ax.set_yticks(range(K))
ax.set_yticklabels([f'Itération {k+1}' for k in range(K)])
ax.set_xticks(np.arange(K) + 0.5)
ax.set_xticklabels([f'Bloc {k+1}' for k in range(K)])
ax.set_xlabel('Blocs de données')
ax.legend(loc='upper right')
ax.set_title(f'Validation croisée à {K} blocs')
ax.set_xlim(0, K)

plt.tight_layout()
<Figure size 1000x300 with 1 Axes>

Le cas particulier K=NK = N (un bloc par exemple) est appelé validation croisée leave-one-out. Elle utilise au maximum les données disponibles, mais est coûteuse en calcul. En pratique, K=5K = 5 ou K=10K = 10 offrent un bon compromis.

Biais inductifs

Il n’existe pas de modèle universel qui fonctionne optimalement pour tous les problèmes. Ce résultat, connu sous le nom de théorème du no free lunch, affirme qu’un algorithme d’apprentissage qui performe bien sur une classe de problèmes performe nécessairement moins bien sur d’autres.

Tout modèle encode des biais inductifs: des hypothèses implicites ou explicites sur la structure du problème. La régression linéaire suppose que la relation entre entrées et sorties est linéaire. Les k plus proches voisins supposent que les points proches dans l’espace des entrées ont des sorties similaires. Les modèles plus complexes, comme les réseaux de neurones, encodent d’autres hypothèses sur la structure des données.

Ces hypothèses sont nécessaires pour que l’apprentissage soit possible. Sans elles, nous n’aurions aucune raison de croire que la performance sur l’échantillon d’entraînement prédit la performance sur de nouvelles données. Le choix du modèle et de ses hypothèses est une décision que l’algorithme ne peut pas prendre seul; elle requiert une connaissance du domaine.

Fonctions de perte de substitution

La perte 0-1 pose un problème pratique. Les méthodes d’optimisation itératives, comme la descente de gradient, requièrent que la fonction objectif soit différentiable. Or la perte 0-1 est constante par morceaux: sa dérivée est nulle presque partout et indéfinie aux points de discontinuité.

Nous contournons ce problème en utilisant des fonctions de perte de substitution: des approximations convexes et différentiables de la perte originale.

Pour la classification binaire, plutôt que de prédire directement une classe, les modèles produisent souvent un score s=f(x)s = f(\mathbf{x}) (un nombre réel). La prédiction de classe se fait ensuite en prenant le signe de ce score: si s>0s > 0, on prédit la classe +1; si s<0s < 0, on prédit la classe -1. La valeur absolue de ss mesure la confiance: plus s|s| est grand, plus le modèle est confiant dans sa prédiction.

Pour la classification binaire avec y{1,+1}y \in \{-1, +1\}, la perte logistique est:

log(y,s)=log(1+eys)\ell_{\text{log}}(y, s) = \log(1 + e^{-y \cdot s})

s=f(x)s = f(\mathbf{x}) est le score produit par le modèle. Cette fonction est convexe et différentiable partout. Lorsque yy et ss ont le même signe (prédiction correcte avec confiance), la perte est faible. Lorsqu’ils ont des signes opposés (erreur), la perte croît linéairement avec l’amplitude de l’erreur.

La perte à charnière (hinge loss) est utilisée dans les machines à vecteurs de support:

hinge(y,s)=max(0,1ys)\ell_{\text{hinge}}(y, s) = \max(0, 1 - y \cdot s)

Cette fonction est convexe mais non différentiable au point ys=1y \cdot s = 1. Elle est nulle lorsque la prédiction est correcte avec une marge suffisante (ys1y \cdot s \geq 1), et croît linéairement sinon.

Ces deux fonctions majorent la perte 0-1: pour tout yy et ss, nous avons 01log\ell_{0-1} \leq \ell_{\text{log}} et 01hinge\ell_{0-1} \leq \ell_{\text{hinge}}. Minimiser ces substituts garantit donc un certain contrôle sur la perte originale.

Source
import numpy as np
import matplotlib.pyplot as plt

# Margin: y * s (positive = correct prediction, negative = error)
margin = np.linspace(-3, 3, 500)

# 0-1 loss: 1 if margin < 0, else 0
loss_01 = (margin < 0).astype(float)

# Logistic loss: log(1 + exp(-margin))
loss_log = np.log(1 + np.exp(-margin))

# Hinge loss: max(0, 1 - margin)
loss_hinge = np.maximum(0, 1 - margin)

fig, ax = plt.subplots(figsize=(8, 5))

ax.plot(margin, loss_01, 'k-', linewidth=2, label='Perte 0-1')
ax.plot(margin, loss_log, 'C0-', linewidth=2, label='Perte logistique')
ax.plot(margin, loss_hinge, 'C1-', linewidth=2, label='Perte à charnière')

ax.axvline(0, color='gray', linestyle=':', alpha=0.5)
ax.axhline(1, color='gray', linestyle=':', alpha=0.3)

ax.set_xlabel(r'Marge $y \cdot s$')
ax.set_ylabel('Perte')
ax.set_xlim(-3, 3)
ax.set_ylim(-0.1, 4)
ax.legend()
ax.set_title('Fonctions de perte de substitution comme bornes supérieures convexes')

# Annotate regions
ax.text(-1.5, 3.5, 'Erreur\n(prédiction incorrecte)', ha='center', fontsize=9, color='gray')
ax.text(1.5, 0.3, 'Correct\n(prédiction juste)', ha='center', fontsize=9, color='gray')

plt.tight_layout()
<Figure size 800x500 with 1 Axes>

La figure montre les trois fonctions de perte en fonction de la marge ysy \cdot s. Une marge positive indique une prédiction correcte (le signe de ss correspond à yy), une marge négative indique une erreur. La perte 0-1 est discontinue au point ys=0y \cdot s = 0. Les pertes logistique et à charnière sont continues et convexes, ce qui permet d’utiliser des méthodes d’optimisation par gradient. Elles majorent partout la perte 0-1.

Le cadre probabiliste

Jusqu’ici, nous avons choisi des fonctions de perte de manière ad hoc: la perte quadratique semble raisonnable pour la régression, la perte logistique pour la classification. Mais d’où viennent ces choix? Existe-t-il un principe unificateur?

Le cadre probabiliste offre une réponse: plutôt que de choisir une perte arbitraire, nous modélisons explicitement comment les données ont été générées. Cette section présente d’abord le cadre général de l’inférence bayésienne, puis développe deux approches concrètes: le maximum de vraisemblance (EMV) et le maximum a posteriori (MAP).

Le cadre bayésien

La statistique bayésienne propose un cadre général pour l’estimation de paramètres. Au lieu d’estimer un point unique, elle caractérise notre incertitude sur les paramètres par une distribution de probabilité.

Le théorème de Bayes nous dit comment mettre à jour nos croyances sur les paramètres θ\boldsymbol{\theta} après avoir observé des données D\mathcal{D}:

p(θD)=p(θ)p(Dθ)p(D)p(\boldsymbol{\theta} | \mathcal{D}) = \frac{p(\boldsymbol{\theta}) \, p(\mathcal{D} | \boldsymbol{\theta})}{p(\mathcal{D})}

Chaque terme a un nom et un rôle précis:

L’a priori encode notre connaissance préalable. Pour une pièce de monnaie, nous pourrions croire que θ\theta est probablement proche de 0.5. L’a posteriori combine cette croyance avec l’évidence des données.

Le prédicteur de Bayes optimal

Commençons par un idéal théorique. Si nous connaissions la vraie distribution conjointe p(x,y)p(\mathbf{x}, y), quelle fonction ff minimiserait le risque?

R(f)=Ep(x,y)[(y,f(x))]=(y,f(x))p(yx)p(x)dydx\mathcal{R}(f) = \mathbb{E}_{p(\mathbf{x}, y)}[\ell(y, f(\mathbf{x}))] = \int \int \ell(y, f(\mathbf{x})) \, p(y | \mathbf{x}) \, p(\mathbf{x}) \, dy \, d\mathbf{x}

Puisque p(x)p(\mathbf{x}) est toujours positif, minimiser cette intégrale revient à minimiser, pour chaque x\mathbf{x}, l’espérance conditionnelle de la perte. Le prédicteur de Bayes optimal est donc:

f(x)=argminy^Ep(yx)[(y,y^)]f^*(\mathbf{x}) = \arg\min_{\hat{y}} \mathbb{E}_{p(y|\mathbf{x})}[\ell(y, \hat{y})]

La réponse dépend de la fonction de perte choisie.

Pour la perte quadratique (y,y^)=(yy^)2\ell(y, \hat{y}) = (y - \hat{y})^2, développons l’espérance:

E[(yy^)2x]=E[y2x]2y^E[yx]+y^2\mathbb{E}[(y - \hat{y})^2 | \mathbf{x}] = \mathbb{E}[y^2 | \mathbf{x}] - 2\hat{y}\mathbb{E}[y | \mathbf{x}] + \hat{y}^2

C’est une fonction quadratique en y^\hat{y}. En dérivant et en posant la dérivée égale à zéro:

y^E[(yy^)2x]=2E[yx]+2y^=0y^=E[yx]\frac{\partial}{\partial \hat{y}} \mathbb{E}[(y - \hat{y})^2 | \mathbf{x}] = -2\mathbb{E}[y | \mathbf{x}] + 2\hat{y} = 0 \quad \Rightarrow \quad \hat{y}^* = \mathbb{E}[y | \mathbf{x}]

Le prédicteur optimal est la moyenne conditionnelle.

Pour la perte 0-1 en classification (y,y^)=1[yy^]\ell(y, \hat{y}) = \mathbf{1}[y \neq \hat{y}], l’espérance est:

E[1[yy^]x]=P(yy^x)=1P(y=y^x)\mathbb{E}[\mathbf{1}[y \neq \hat{y}] | \mathbf{x}] = P(y \neq \hat{y} | \mathbf{x}) = 1 - P(y = \hat{y} | \mathbf{x})

Minimiser cette quantité revient à maximiser P(y=y^x)P(y = \hat{y} | \mathbf{x}), donc à choisir la classe la plus probable:

y^=argmaxcp(y=cx)\hat{y}^* = \arg\max_c \, p(y = c | \mathbf{x})

Le prédicteur optimal est le mode conditionnel. Chaque perte définit son propre prédicteur optimal.

Ce prédicteur est un repère théorique: aucun algorithme ne peut faire mieux, car il suppose l’accès à la vraie distribution. La différence entre le risque d’un prédicteur appris et ce risque de Bayes mesure ce que nous perdons en ne connaissant pas la vraie distribution.

Prédiction bayésienne et distribution prédictive a posteriori

En pratique, nous ne connaissons pas p(yx)p(y|\mathbf{x}). Nous avons un modèle paramétrique p(yx,θ)p(y|\mathbf{x}, \boldsymbol{\theta}) et une distribution a posteriori p(θD)p(\boldsymbol{\theta}|\mathcal{D}) sur les paramètres. L’approche bayésienne complète consiste à moyenner les prédictions sur tous les paramètres possibles, pondérés par leur probabilité a posteriori:

p(yx,D)=p(yx,θ)p(θD)dθp(y|\mathbf{x}, \mathcal{D}) = \int p(y|\mathbf{x}, \boldsymbol{\theta}) \, p(\boldsymbol{\theta}|\mathcal{D}) \, d\boldsymbol{\theta}

Cette distribution prédictive a posteriori intègre l’incertitude sur les paramètres. Elle ne s’engage pas sur une valeur unique de θ\boldsymbol{\theta}, mais considère toutes les valeurs plausibles.

Le problème: cette intégrale est rarement calculable analytiquement. Elle nécessite de sommer sur un espace de paramètres de grande dimension, ce qui est coûteux ou impossible en pratique. C’est pourquoi nous recourons souvent à des estimateurs ponctuels: plutôt que d’intégrer sur tous les θ\boldsymbol{\theta}, nous en choisissons un seul, comme l’EMV ou le MAP.

Utilité du modèle probabiliste

Si nous finissons souvent par utiliser un estimateur ponctuel, pourquoi adopter le cadre probabiliste? Plusieurs raisons:

  1. Justifier la fonction de perte: La perte quadratique découle naturellement de l’hypothèse de bruit gaussien. La perte logarithmique vient du principe de maximum de vraisemblance. Le cadre probabiliste explique pourquoi ces choix sont raisonnables.

  2. Quantifier l’incertitude: Au-delà de la prédiction ponctuelle y^=f(x;θ^)\hat{y} = f(\mathbf{x}; \hat{\boldsymbol{\theta}}), nous pouvons donner un intervalle de prédiction. Sous un modèle gaussien, yy a environ 95% de chances de tomber dans [f(x)2σ,f(x)+2σ][f(\mathbf{x}) - 2\sigma, f(\mathbf{x}) + 2\sigma].

  3. Comparer des modèles: La vraisemblance marginale p(D)p(\mathcal{D}) permet de comparer des modèles de complexités différentes, pénalisant automatiquement les modèles trop complexes.

  4. Ouvrir la porte à l’inférence complète: Quand les ressources le permettent (méthodes de Monte Carlo, inférence variationnelle), nous pouvons approximer la distribution prédictive complète plutôt que de nous limiter à un point.

Maximum de vraisemblance

Le maximum de vraisemblance est la première approche concrète dans le cadre probabiliste: nous cherchons les paramètres qui rendent nos observations les plus probables.

Construction de la vraisemblance

Supposons que nous avons un modèle paramétrique p(yx;θ)p(y|\mathbf{x}; \boldsymbol{\theta}) qui, pour chaque entrée x\mathbf{x} et choix de paramètres θ\boldsymbol{\theta}, définit une distribution sur les sorties possibles yy. Par exemple, en régression, ce pourrait être une gaussienne centrée sur f(x;θ)f(\mathbf{x}; \boldsymbol{\theta}).

Considérons un seul exemple (x1,y1)(\mathbf{x}_1, y_1). Pour des paramètres θ\boldsymbol{\theta} fixés, nous pouvons évaluer p(y1x1;θ)p(y_1 | \mathbf{x}_1; \boldsymbol{\theta}): la probabilité (ou densité) que le modèle assigne à l’observation y1y_1. Si cette valeur est élevée, les paramètres θ\boldsymbol{\theta} “expliquent bien” cette observation. Si elle est faible, y1y_1 est une valeur improbable sous ce modèle.

Avec deux exemples indépendants (x1,y1)(\mathbf{x}_1, y_1) et (x2,y2)(\mathbf{x}_2, y_2), la probabilité conjointe est le produit:

p(y1,y2x1,x2;θ)=p(y1x1;θ)p(y2x2;θ)p(y_1, y_2 | \mathbf{x}_1, \mathbf{x}_2; \boldsymbol{\theta}) = p(y_1 | \mathbf{x}_1; \boldsymbol{\theta}) \cdot p(y_2 | \mathbf{x}_2; \boldsymbol{\theta})

Avec NN exemples indépendants, nous obtenons la vraisemblance:

L(θ)=i=1Np(yixi;θ)\mathcal{L}(\boldsymbol{\theta}) = \prod_{i=1}^N p(y_i | \mathbf{x}_i; \boldsymbol{\theta})

Cette quantité est une fonction de θ\boldsymbol{\theta}. Elle répond à la question: pour ce choix de paramètres, quelle est la probabilité d’avoir observé exactement ces données?

Pourquoi maximiser?

Si L(θA)>L(θB)\mathcal{L}(\boldsymbol{\theta}_A) > \mathcal{L}(\boldsymbol{\theta}_B), alors les données observées sont plus probables sous θA\boldsymbol{\theta}_A que sous θB\boldsymbol{\theta}_B. Les paramètres θA\boldsymbol{\theta}_A rendent les observations moins “surprenantes”.

L’estimateur du maximum de vraisemblance (EMV, ou MLE pour maximum likelihood estimator en anglais) choisit les paramètres qui maximisent cette probabilité:

θ^EMV=argmaxθL(θ)=argmaxθi=1Np(yixi;θ)\hat{\boldsymbol{\theta}}_{\text{EMV}} = \arg\max_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) = \arg\max_{\boldsymbol{\theta}} \prod_{i=1}^N p(y_i | \mathbf{x}_i; \boldsymbol{\theta})

C’est le choix de paramètres sous lequel nos données sont les plus “attendues”.

Du produit à la somme

En pratique, multiplier NN probabilités (souvent petites) pose des problèmes numériques: le résultat devient rapidement trop petit pour être représenté par un ordinateur. Le logarithme résout ce problème: il transforme le produit en somme et, comme c’est une fonction croissante, il ne change pas le maximiseur:

logL(θ)=i=1Nlogp(yixi;θ)\log \mathcal{L}(\boldsymbol{\theta}) = \sum_{i=1}^N \log p(y_i | \mathbf{x}_i; \boldsymbol{\theta})

Pour l’optimisation, nous préférons minimiser plutôt que maximiser (par convention). La log-vraisemblance négative (NLV, ou NLL pour negative log-likelihood en anglais) est notre fonction objectif:

NLV(θ)=i=1Nlogp(yixi;θ)\text{NLV}(\boldsymbol{\theta}) = -\sum_{i=1}^N \log p(y_i | \mathbf{x}_i; \boldsymbol{\theta})

Remarquez la structure: c’est une somme sur les exemples d’une quantité logp(yixi;θ)-\log p(y_i | \mathbf{x}_i; \boldsymbol{\theta}) qui dépend de chaque observation. Cette quantité joue le rôle d’une fonction de perte. Le maximum de vraisemblance est donc un cas particulier de la minimisation du risque empirique, où la perte est définie par le modèle probabiliste lui-même.

Régression avec bruit gaussien: d’où vient la perte quadratique?

Appliquons ce principe à la régression. Le modèle de génération des données suppose que la sortie observée est la prédiction “vraie” du modèle, corrompue par un bruit aléatoire gaussien:

y=f(x;θ)+ε,εN(0,σ2)y = f(\mathbf{x}; \boldsymbol{\theta}) + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma^2)

Ce modèle dit que si nous connaissions les vrais paramètres θ\boldsymbol{\theta} et que nous mesurions yy pour un x\mathbf{x} donné, nous obtiendrions f(x;θ)f(\mathbf{x}; \boldsymbol{\theta}) plus ou moins σ\sigma la plupart du temps.

La distribution conditionnelle qui en découle est:

p(yx;θ)=12πσ2exp((yf(x;θ))22σ2)p(y|\mathbf{x}; \boldsymbol{\theta}) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(y - f(\mathbf{x}; \boldsymbol{\theta}))^2}{2\sigma^2}\right)

Calculons la log-vraisemblance négative:

NLV(θ)=i=1Nlogp(yixi;θ)=12σ2i=1N(yif(xi;θ))2+N2log(2πσ2)\text{NLV}(\boldsymbol{\theta}) = -\sum_{i=1}^N \log p(y_i | \mathbf{x}_i; \boldsymbol{\theta}) = \frac{1}{2\sigma^2} \sum_{i=1}^N (y_i - f(\mathbf{x}_i; \boldsymbol{\theta}))^2 + \frac{N}{2}\log(2\pi\sigma^2)

Le second terme ne dépend pas de θ\boldsymbol{\theta}. Minimiser la NLV revient donc exactement à minimiser la somme des erreurs quadratiques.

C’est un résultat fondamental: la perte quadratique n’est pas un choix arbitraire. Elle découle naturellement de l’hypothèse que les erreurs de mesure suivent une loi gaussienne. Le maximum de vraisemblance sous bruit gaussien coïncide avec les moindres carrés.

Dans ce modèle, nous avons supposé que la variance σ2\sigma^2 est constante pour toutes les entrées x\mathbf{x}. C’est ce qu’on appelle la régression homoscédastique (du grec homos, même, et skedasis, dispersion). C’est l’hypothèse standard en régression linéaire.

En pratique, l’incertitude peut varier selon l’entrée. Par exemple, les mesures à haute vitesse peuvent être plus bruitées que celles à basse vitesse. La régression hétéroscédastique modélise cette variation en faisant dépendre la variance de x\mathbf{x}:

p(yx;θ)=N(yfμ(x;θ),fσ(x;θ)2)p(y|\mathbf{x}; \boldsymbol{\theta}) = \mathcal{N}(y | f_\mu(\mathbf{x}; \boldsymbol{\theta}), f_\sigma(\mathbf{x}; \boldsymbol{\theta})^2)

fμf_\mu prédit la moyenne et fσf_\sigma prédit l’écart-type. Ce modèle est plus flexible mais requiert d’apprendre des paramètres supplémentaires.

Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.stats import norm
from IPython.display import HTML

# Générer des données synthétiques
np.random.seed(42)
N = 100
x_data = np.random.uniform(0.5, 9.5, N)
f_mu = lambda x: 0.5 * x + 1

# Homoscédastique: variance constante
sigma_homo = 0.7
y_homo = f_mu(x_data) + np.random.normal(0, sigma_homo, N)

# Hétéroscédastique: variance croissante
f_sigma = lambda x: 0.3 + 0.12 * x
y_hetero = f_mu(x_data) + np.random.normal(0, f_sigma(x_data))

# Configuration de la figure
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
x_line = np.linspace(0, 10, 100)
y_pdf_range = np.linspace(-3, 9, 200)
scale = 2.5  # échelle pour afficher les PDFs

def init():
    for ax in axes:
        ax.clear()
    return []

def animate(frame):
    x_current = 0.5 + frame * 9 / 59  # balayer de 0.5 à 9.5
    
    for idx, (ax, y_data, title, color, get_sigma) in enumerate([
        (axes[0], y_homo, 'Régression homoscédastique', 'steelblue', lambda x: sigma_homo),
        (axes[1], y_hetero, 'Régression hétéroscédastique', 'coral', f_sigma)
    ]):
        ax.clear()
        
        # Données et ligne de régression
        ax.scatter(x_data, y_data, alpha=0.4, s=20, c='gray', zorder=1)
        ax.plot(x_line, f_mu(x_line), 'k-', linewidth=2, zorder=2)
        
        # Gaussienne à la position actuelle
        mu = f_mu(x_current)
        sigma = get_sigma(x_current)
        pdf = norm.pdf(y_pdf_range, mu, sigma)
        
        # Afficher la gaussienne "horizontalement"
        ax.fill_betweenx(y_pdf_range, x_current, x_current + scale * pdf, 
                         alpha=0.5, color=color, zorder=3)
        ax.plot(x_current + scale * pdf, y_pdf_range, color=color, linewidth=2, zorder=4)
        
        # Ligne verticale indiquant la position
        ax.axvline(x_current, color=color, linestyle='--', alpha=0.5, linewidth=1)
        
        # Point sur la courbe de régression
        ax.scatter([x_current], [mu], color='black', s=50, zorder=5)
        
        # Bande ±2σ
        ax.fill_between([x_current - 0.1, x_current + 0.1], 
                        [mu - 2*sigma, mu - 2*sigma], 
                        [mu + 2*sigma, mu + 2*sigma],
                        alpha=0.2, color=color, zorder=0)
        
        ax.set_xlim(-0.5, 12)
        ax.set_ylim(-2, 8)
        ax.set_xlabel(r'$x$', fontsize=11)
        ax.set_ylabel(r'$y$', fontsize=11)
        sigma_label = r'$\sigma^2$ constant' if idx == 0 else r'$\sigma^2(x)$ variable'
        ax.set_title(f'{title}\n{sigma_label}', fontsize=11)
    
    fig.tight_layout()
    return []

anim = FuncAnimation(fig, animate, init_func=init, frames=60, interval=80, blit=True)
anim.save('_static/regression_scedasticity.gif', writer='pillow', fps=12, dpi=100)
plt.close()

# Afficher le GIF
from IPython.display import Image
Image(filename='_static/regression_scedasticity.gif')
<IPython.core.display.Image object>

L’animation illustre la différence fondamentale entre les deux modèles. À chaque position xx, la distribution conditionnelle p(yx)p(y|x) est une gaussienne (la “cloche” colorée) centrée sur la courbe de régression fμ(x)f_\mu(x). Dans le cas homoscédastique (gauche), la cloche garde la même largeur partout. Dans le cas hétéroscédastique (droite), la largeur varie avec xx. Ici, l’incertitude augmente vers la droite, ce qui se traduit par une dispersion plus grande des points.

Régression linéaire homoscédastique

Appliquons maintenant le maximum de vraisemblance au cas le plus courant: la régression linéaire homoscédastique. Le modèle probabiliste est:

p(yx;θ,σ2)=N(yθx,σ2)p(y | \mathbf{x}; \boldsymbol{\theta}, \sigma^2) = \mathcal{N}(y | \boldsymbol{\theta}^\top \mathbf{x}, \sigma^2)

La fonction de moyenne est linéaire: fμ(x;θ)=θx=j=0dθjxjf_\mu(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \mathbf{x} = \sum_{j=0}^d \theta_j x_j, où nous avons absorbé le biais en posant x0=1x_0 = 1. Le vecteur θRd+1\boldsymbol{\theta} \in \mathbb{R}^{d+1} contient donc le biais θ0\theta_0 et les poids θ1,,θd\theta_1, \ldots, \theta_d. La variance σ2\sigma^2 est constante.

Formulation matricielle

Avec NN observations {(xi,yi)}i=1N\{(\mathbf{x}_i, y_i)\}_{i=1}^N, nous pouvons écrire le modèle sous forme matricielle. Définissons la matrice de conception (design matrix) XRN×(d+1)\mathbf{X} \in \mathbb{R}^{N \times (d+1)} dont chaque ligne contient une observation augmentée d’un 1 pour le biais:

X=(1x11x12x1d1x21x22x2d1xN1xN2xNd),y=(y1y2yN),θ=(θ0θ1θd)\mathbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1d} \\ 1 & x_{21} & x_{22} & \cdots & x_{2d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & x_{N2} & \cdots & x_{Nd} \end{pmatrix}, \quad \mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix}, \quad \boldsymbol{\theta} = \begin{pmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_d \end{pmatrix}

Le vecteur des prédictions est y^=Xθ\hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\theta}, et le vecteur des résidus est r=yXθ\mathbf{r} = \mathbf{y} - \mathbf{X}\boldsymbol{\theta}.

La somme des carrés des résidus (residual sum of squares, RSS) s’écrit:

RSS(θ)=i=1N(yiθxi)2=yXθ22=(yXθ)(yXθ)\text{RSS}(\boldsymbol{\theta}) = \sum_{i=1}^N (y_i - \boldsymbol{\theta}^\top \mathbf{x}_i)^2 = \|\mathbf{y} - \mathbf{X}\boldsymbol{\theta}\|_2^2 = (\mathbf{y} - \mathbf{X}\boldsymbol{\theta})^\top (\mathbf{y} - \mathbf{X}\boldsymbol{\theta})

Comme nous l’avons vu, minimiser la NLV sous bruit gaussien homoscédastique revient à minimiser le RSS (ou de manière équivalente, l’EQM = RSS/NN).

Solution analytique: les équations normales

Pour trouver le minimum, nous calculons le gradient du RSS par rapport à θ\boldsymbol{\theta} et l’égalons à zéro. En développant:

RSS(θ)=yy2θXy+θXXθ\text{RSS}(\boldsymbol{\theta}) = \mathbf{y}^\top \mathbf{y} - 2\boldsymbol{\theta}^\top \mathbf{X}^\top \mathbf{y} + \boldsymbol{\theta}^\top \mathbf{X}^\top \mathbf{X} \boldsymbol{\theta}

Le gradient est:

θRSS(θ)=2Xy+2XXθ\nabla_{\boldsymbol{\theta}} \text{RSS}(\boldsymbol{\theta}) = -2\mathbf{X}^\top \mathbf{y} + 2\mathbf{X}^\top \mathbf{X} \boldsymbol{\theta}

En posant ce gradient égal à zéro, nous obtenons les équations normales:

XXθ=Xy\mathbf{X}^\top \mathbf{X} \boldsymbol{\theta} = \mathbf{X}^\top \mathbf{y}

Si la matrice XX\mathbf{X}^\top \mathbf{X} est inversible, la solution unique est:

θ^EMV=(XX)1Xy\hat{\boldsymbol{\theta}}_{\text{EMV}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}

Cette solution porte le nom d’estimateur des moindres carrés ordinaires (MCO, ou ordinary least squares, OLS). Elle est exactement équivalente à l’EMV sous l’hypothèse de bruit gaussien homoscédastique.

Conditions d’existence de la solution

La matrice XX\mathbf{X}^\top \mathbf{X} est de taille (d+1)×(d+1)(d+1) \times (d+1). Elle est inversible si et seulement si X\mathbf{X} est de rang plein colonne, c’est-à-dire si les colonnes de X\mathbf{X} sont linéairement indépendantes. Cela requiert:

  1. Plus d’observations que de paramètres: Nd+1N \geq d + 1

  2. Pas de colinéarité parfaite: aucune caractéristique ne doit être une combinaison linéaire exacte des autres

Quand XX\mathbf{X}^\top \mathbf{X} est mal conditionnée (presque singulière), de petites perturbations dans les données peuvent causer de grandes variations dans la solution. C’est l’un des problèmes que la régularisation permet de résoudre.

Source
import numpy as np
import matplotlib.pyplot as plt

# Données synthétiques: relation linéaire avec bruit
np.random.seed(42)
N = 30
x = np.random.uniform(0, 10, N)
y_true = 2.5 * x + 3  # vraie relation: y = 2.5x + 3
y = y_true + np.random.normal(0, 3, N)  # ajout de bruit gaussien

# Construction de la matrice de conception (avec colonne de 1 pour le biais)
X = np.column_stack([np.ones(N), x])

# Solution OLS: w = (X^T X)^{-1} X^T y
XtX = X.T @ X
Xty = X.T @ y
w_ols = np.linalg.solve(XtX, Xty)
b_ols, slope_ols = w_ols[0], w_ols[1]

# Prédictions
x_grid = np.linspace(0, 10, 100)
y_pred = b_ols + slope_ols * x_grid

# Visualisation
fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Gauche: données et droite ajustée
ax = axes[0]
ax.scatter(x, y, alpha=0.7, label='Observations')
ax.plot(x_grid, y_pred, 'k-', linewidth=2, label=f'MCO: $y = {slope_ols:.2f}x + {b_ols:.2f}$')
ax.plot(x_grid, 2.5 * x_grid + 3, 'g--', alpha=0.5, label='Vraie relation')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.legend()
ax.set_title('Régression linéaire par moindres carrés')

# Droite: résidus
ax = axes[1]
y_fitted = b_ols + slope_ols * x
residuals = y - y_fitted
ax.stem(x, residuals, basefmt=' ', linefmt='C0-', markerfmt='C0o')
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel('$x$')
ax.set_ylabel(r'Résidu $r_i = y_i - \hat{y}_i$')
ax.set_title(f'Résidus (RSS = {np.sum(residuals**2):.1f})')

plt.tight_layout()
<Figure size 1000x400 with 2 Axes>

La figure de gauche montre les données et la droite ajustée par moindres carrés. La figure de droite montre les résidus: les écarts entre les observations et les prédictions. L’EMV minimise la somme des carrés de ces résidus.

Exemple: pharmacocinétique

L’EMV s’applique à des modèles non linéaires. Considérons la concentration d’un médicament dans le sang après administration orale. Les données suivantes proviennent d’une étude sur la théophylline, un bronchodilatateur:

Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# Données pharmacocinétiques: théophylline, sujet 1 (Boeckmann et al., 1994)
time = np.array([0, 0.25, 0.57, 1.12, 2.02, 3.82, 5.10, 7.03, 9.05, 12.12, 24.37])
conc = np.array([0.74, 2.84, 6.57, 10.50, 9.66, 8.58, 8.36, 7.47, 6.89, 5.94, 3.28])

# Model: C(t) = C0 * exp(-k * t) for t > t_peak
# We'll fit on the decay phase (after peak)
peak_idx = np.argmax(conc)
t_decay = time[peak_idx:]
c_decay = conc[peak_idx:]

# MLE: minimize NLL under Gaussian noise
def neg_log_likelihood(params, t, c):
    C0, k, sigma = params
    if sigma <= 0 or k <= 0 or C0 <= 0:
        return np.inf
    pred = C0 * np.exp(-k * (t - t[0]))
    nll = 0.5 * len(t) * np.log(2 * np.pi * sigma**2)
    nll += 0.5 * np.sum((c - pred)**2) / sigma**2
    return nll

# Initial guess and optimization
x0 = [c_decay[0], 0.1, 1.0]
result = minimize(neg_log_likelihood, x0, args=(t_decay, c_decay), method='Nelder-Mead')
C0_mle, k_mle, sigma_mle = result.x

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Left: data and fit
ax = axes[0]
ax.scatter(time, conc, s=50, zorder=5, label='Observations')
t_grid = np.linspace(time[peak_idx], 25, 100)
ax.plot(t_grid, C0_mle * np.exp(-k_mle * (t_grid - time[peak_idx])), 'k--', 
        label=f'EMV: $C_0$={C0_mle:.1f}, $k$={k_mle:.2f}')
ax.axvline(time[peak_idx], color='gray', linestyle=':', alpha=0.5)
ax.set_xlabel('Temps (h)')
ax.set_ylabel('Concentration (mg/L)')
ax.legend()
ax.set_title('Concentration plasmatique de théophylline')

# Right: residuals
ax = axes[1]
pred_decay = C0_mle * np.exp(-k_mle * (t_decay - t_decay[0]))
residuals = c_decay - pred_decay
ax.stem(t_decay, residuals, basefmt=' ')
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel('Temps (h)')
ax.set_ylabel('Résidu (mg/L)')
ax.set_title(rf'$\sigma$ estimé: {sigma_mle:.2f} mg/L')

plt.tight_layout()
<Figure size 1000x400 with 2 Axes>

Le modèle C(t)=C0ektC(t) = C_0 e^{-kt} décrit la décroissance exponentielle après le pic de concentration. Les paramètres C0C_0 (concentration initiale) et kk (constante d’élimination) sont estimés par maximum de vraisemblance sous l’hypothèse d’un bruit gaussien. Cette approche est identique à celle des moindres carrés, mais elle fournit également une estimation de l’écart-type du bruit σ\sigma.

Classification binaire

La perte 0-1 pour la classification est discontinue, ce qui empêche l’utilisation de méthodes de gradient. La fonction sigmoïde σ(z)=1/(1+ez)\sigma(z) = 1/(1 + e^{-z}) contourne ce problème: c’est une approximation lisse de la fonction échelon (step function). Elle transforme n’importe quel score réel en une valeur dans l’intervalle (0,1)(0, 1), que nous pouvons interpréter comme une probabilité.

Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Create figure
fig, ax = plt.subplots(figsize=(8, 5))

# Define x range
z = np.linspace(-4, 4, 200)

# Step function (Heaviside)
step = (z >= 0).astype(float)

# Sigmoid function with temperature parameter
def sigmoid(z, alpha=1):
    return 1 / (1 + np.exp(-alpha * z))

# Initialize plot
line_step, = ax.plot(z, step, 'k--', linewidth=2, label='Fonction échelon', alpha=0.7)
line_sigmoid, = ax.plot([], [], 'b-', linewidth=2, label='Sigmoïde $\\sigma(\\alpha z)$')
ax.axhline(0.5, color='gray', linestyle=':', alpha=0.5, linewidth=1)
ax.axvline(0, color='gray', linestyle=':', alpha=0.5, linewidth=1)
ax.set_xlim(-4, 4)
ax.set_ylim(-0.1, 1.1)
ax.set_xlabel('$z$')
ax.set_ylabel('$\\sigma(\\alpha z)$')
ax.set_title('Approximation de la fonction échelon par la sigmoïde')
ax.legend(loc='best')
ax.grid(True, alpha=0.3)

# Animation function
def animate(frame):
    # Alpha increases from 0.5 to 10
    alpha = 0.5 + (frame / 100) * 9.5
    y = sigmoid(z, alpha)
    line_sigmoid.set_data(z, y)
    ax.set_title(f'Approximation de la fonction échelon par la sigmoïde ($\\alpha = {alpha:.2f}$)')
    return line_sigmoid,

# Create animation
anim = FuncAnimation(fig, animate, frames=100, interval=50, blit=True, repeat=True)
anim.save('_static/sigmoid_approximation.gif', writer='pillow', fps=20, dpi=100)
plt.close()

# Afficher le GIF
from IPython.display import Image
Image(filename='_static/sigmoid_approximation.gif')
<IPython.core.display.Image object>

L’animation montre comment la sigmoïde σ(αz)\sigma(\alpha z) se rapproche de la fonction échelon lorsque le paramètre α\alpha augmente. Pour α=1\alpha = 1, la sigmoïde est douce; pour α\alpha grand, elle devient presque aussi abrupte que la fonction échelon, tout en restant différentiable.

Cette interprétation probabiliste n’est pas qu’une astuce numérique. Elle correspond exactement à modéliser YXY | \mathbf{X} par une distribution de Bernoulli dont le paramètre dépend de l’entrée.

Pour la classification binaire avec y{0,1}y \in \{0, 1\}, nous modélisons la probabilité de la classe positive par:

p(y=1x;θ)=σ(f(x;θ))=11+ef(x;θ)p(y = 1 | \mathbf{x}; \boldsymbol{\theta}) = \sigma(f(\mathbf{x}; \boldsymbol{\theta})) = \frac{1}{1 + e^{-f(\mathbf{x}; \boldsymbol{\theta})}}

σ\sigma est la fonction sigmoïde et f(x;θ)f(\mathbf{x}; \boldsymbol{\theta}) est le logit (ou log-odds), le score brut du modèle avant transformation. Le logit est le logarithme du rapport des probabilités: logp(y=1x)p(y=0x)=logp1p\log \frac{p(y=1|\mathbf{x})}{p(y=0|\mathbf{x})} = \log \frac{p}{1-p}. La distribution conditionnelle suit une loi de Bernoulli:

p(yx;θ)=σ(f(x;θ))y(1σ(f(x;θ)))1yp(y|\mathbf{x}; \boldsymbol{\theta}) = \sigma(f(\mathbf{x}; \boldsymbol{\theta}))^y (1 - \sigma(f(\mathbf{x}; \boldsymbol{\theta})))^{1-y}

La log-vraisemblance négative est:

NLV(θ)=i=1N[yilogσ(f(xi;θ))+(1yi)log(1σ(f(xi;θ)))]\text{NLV}(\boldsymbol{\theta}) = -\sum_{i=1}^N \left[ y_i \log \sigma(f(\mathbf{x}_i; \boldsymbol{\theta})) + (1-y_i) \log(1 - \sigma(f(\mathbf{x}_i; \boldsymbol{\theta}))) \right]

Cette quantité est l’entropie croisée binaire. Elle correspond à la perte logistique, à une reparamétrisation près.

Classification multiclasse

Pour la classification avec CC classes (C>2C > 2), nous généralisons le modèle binaire en utilisant la distribution catégorielle (ou multinomiale). Au lieu de modéliser une seule probabilité p(y=1x)p(y=1|\mathbf{x}), nous modélisons un vecteur de probabilités π(x)=[π1(x),,πC(x)]\boldsymbol{\pi}(\mathbf{x}) = [\pi_1(\mathbf{x}), \ldots, \pi_C(\mathbf{x})]πc(x)=p(y=cx)\pi_c(\mathbf{x}) = p(y=c|\mathbf{x}) et c=1Cπc(x)=1\sum_{c=1}^C \pi_c(\mathbf{x}) = 1.

Pour transformer les scores bruts du modèle en probabilités, nous utilisons la fonction softmax:

πc(x;θ)=exp(fc(x;θ))j=1Cexp(fj(x;θ))\pi_c(\mathbf{x}; \boldsymbol{\theta}) = \frac{\exp(f_c(\mathbf{x}; \boldsymbol{\theta}))}{\sum_{j=1}^C \exp(f_j(\mathbf{x}; \boldsymbol{\theta}))}

fc(x;θ)f_c(\mathbf{x}; \boldsymbol{\theta}) est le score pour la classe cc. La fonction softmax généralise la sigmoïde au cas multiclasse: elle transforme CC scores réels en un vecteur de probabilités qui somme à 1.

La distribution conditionnelle suit une loi catégorielle:

p(yx;θ)=c=1Cπc(x;θ)1[y=c]p(y|\mathbf{x}; \boldsymbol{\theta}) = \prod_{c=1}^C \pi_c(\mathbf{x}; \boldsymbol{\theta})^{\mathbf{1}[y = c]}

1[y=c]\mathbf{1}[y = c] vaut 1 si y=cy = c et 0 sinon. En utilisant l’encodage one-hot y=[1[y=1],,1[y=C]]\mathbf{y} = [\mathbf{1}[y=1], \ldots, \mathbf{1}[y=C]]^\top, cette expression devient:

p(yx;θ)=c=1Cπc(x;θ)ycp(y|\mathbf{x}; \boldsymbol{\theta}) = \prod_{c=1}^C \pi_c(\mathbf{x}; \boldsymbol{\theta})^{y_c}

La log-vraisemblance négative est:

NLV(θ)=i=1Nc=1Cyiclogπc(xi;θ)\text{NLV}(\boldsymbol{\theta}) = -\sum_{i=1}^N \sum_{c=1}^C y_{ic} \log \pi_c(\mathbf{x}_i; \boldsymbol{\theta})

yic=1[yi=c]y_{ic} = \mathbf{1}[y_i = c]. Cette quantité est l’entropie croisée multiclasse. Elle généralise l’entropie croisée binaire au cas où il y a plus de deux classes.

Pour la classification binaire avec C=2C=2, le softmax se réduit à la sigmoïde. En effet, si nous définissons s=f1(x)f2(x)s = f_1(\mathbf{x}) - f_2(\mathbf{x}), alors:

π1=ef1ef1+ef2=11+e(f1f2)=σ(s)\pi_1 = \frac{e^{f_1}}{e^{f_1} + e^{f_2}} = \frac{1}{1 + e^{-(f_1 - f_2)}} = \sigma(s)

Le modèle binaire et le modèle multiclasse partagent donc la même structure probabiliste, avec la distribution catégorielle comme généralisation naturelle de la distribution de Bernoulli.

Maximum a posteriori

Plutôt que de travailler avec la distribution a posteriori complète (ce qui peut être coûteux), nous pouvons chercher son mode: la valeur des paramètres la plus probable a posteriori. C’est l’estimateur du maximum a posteriori (MAP):

θ^MAP=argmaxθp(θD)=argmaxθp(θ)p(Dθ)\hat{\boldsymbol{\theta}}_{\text{MAP}} = \arg\max_{\boldsymbol{\theta}} p(\boldsymbol{\theta} | \mathcal{D}) = \arg\max_{\boldsymbol{\theta}} p(\boldsymbol{\theta}) \, p(\mathcal{D} | \boldsymbol{\theta})

Le dénominateur p(D)p(\mathcal{D}) ne dépend pas de θ\boldsymbol{\theta} et peut être ignoré pour l’optimisation. En passant au logarithme:

θ^MAP=argmaxθ[logp(Dθ)+logp(θ)]\hat{\boldsymbol{\theta}}_{\text{MAP}} = \arg\max_{\boldsymbol{\theta}} \left[ \log p(\mathcal{D} | \boldsymbol{\theta}) + \log p(\boldsymbol{\theta}) \right]

Cette expression révèle une structure familière. Si nous posons C(θ)=logp(θ)C(\boldsymbol{\theta}) = -\log p(\boldsymbol{\theta}), nous obtenons:

θ^MAP=argminθ[NLV(θ)+C(θ)]\hat{\boldsymbol{\theta}}_{\text{MAP}} = \arg\min_{\boldsymbol{\theta}} \left[ \text{NLV}(\boldsymbol{\theta}) + C(\boldsymbol{\theta}) \right]

C’est exactement la forme du risque empirique régularisé. La régularisation correspond à l’ajout d’un a priori sur les paramètres. Le terme de régularisation C(θ)C(\boldsymbol{\theta}) est le logarithme négatif de la distribution a priori.

Le maximum de vraisemblance comme cas particulier

Que se passe-t-il si nous n’avons aucune préférence a priori sur les paramètres? Cela correspond à un a priori uniforme (ou constant): p(θ)=constantep(\boldsymbol{\theta}) = \text{constante}.

Dans ce cas, logp(θ)\log p(\boldsymbol{\theta}) est une constante qui n’affecte pas l’optimisation, et le MAP se réduit à l’EMV:

θ^MAP=θ^EMVquand p(θ)=constante\hat{\boldsymbol{\theta}}_{\text{MAP}} = \hat{\boldsymbol{\theta}}_{\text{EMV}} \quad \text{quand } p(\boldsymbol{\theta}) = \text{constante}

L’EMV est donc un cas particulier du MAP: celui où nous supposons implicitement que toutes les valeurs de paramètres sont également plausibles avant d’observer les données. Cette perspective unifie les deux approches dans un même cadre.

Limites de l’a priori uniforme

L’a priori uniforme (et donc l’EMV) peut être problématique quand les données sont peu nombreuses. Considérons l’estimation de la probabilité θ\theta qu’une pièce tombe sur face.

Supposons que nous lancions la pièce 3 fois et obtenions 3 faces. L’estimateur du maximum de vraisemblance pour une distribution de Bernoulli est:

θ^EMV=N1N0+N1=30+3=1\hat{\theta}_{\text{EMV}} = \frac{N_1}{N_0 + N_1} = \frac{3}{0 + 3} = 1

N1N_1 est le nombre de faces et N0N_0 le nombre de piles. Cette estimation dit que la probabilité d’obtenir face est de 100%. Si nous utilisions ce modèle pour prédire de futurs lancers, nous prédirons toujours face, ce qui est peu plausible pour une vraie pièce.

Le problème est que l’EMV (avec son a priori uniforme implicite) dispose de suffisamment de flexibilité pour reproduire parfaitement les données d’entraînement, même quand celles-ci sont peu nombreuses ou non représentatives. Un a priori informatif peut atténuer ce problème.

Source
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# MLE vs MAP for Bernoulli with few observations
fig, axes = plt.subplots(1, 3, figsize=(12, 4))

# Different sample sizes
samples_list = [
    [1, 1, 1],           # 3 heads
    [1, 1, 1, 0],        # 3 heads, 1 tail
    [1, 1, 1, 0, 0, 1, 1, 0, 1, 1]  # 7 heads, 3 tails
]

theta_grid = np.linspace(0.001, 0.999, 200)

for ax, samples in zip(axes, samples_list):
    n1 = sum(samples)  # heads
    n0 = len(samples) - n1  # tails
    
    # MLE
    theta_mle = n1 / (n0 + n1)
    
    # Likelihood (unnormalized)
    likelihood = theta_grid**n1 * (1 - theta_grid)**n0
    likelihood = likelihood / likelihood.max()
    
    ax.plot(theta_grid, likelihood, 'b-', linewidth=2, label='Vraisemblance')
    ax.axvline(theta_mle, color='b', linestyle='--', alpha=0.7,
               label=f'EMV: {theta_mle:.2f}')
    ax.axvline(0.5, color='gray', linestyle=':', alpha=0.5, label=r'$\theta = 0.5$')
    
    ax.set_xlabel(r'$\theta$')
    ax.set_ylabel('Vraisemblance (normalisée)')
    ax.set_title(f'{n1} faces, {n0} piles (N={len(samples)})')
    ax.legend(fontsize=8)
    ax.set_xlim(0, 1)

plt.tight_layout()
<Figure size 1200x400 with 3 Axes>

La figure montre la vraisemblance pour différents échantillons. Avec seulement 3 observations (toutes faces), la vraisemblance est maximale à θ=1\theta = 1. En augmentant la taille de l’échantillon, l’estimation devient plus raisonnable. Voyons comment un a priori non uniforme peut aider.

Exemple: lissage de Laplace

Revenons à notre exemple de la pièce de monnaie. Utilisons un a priori Beta sur θ\theta:

p(θ)=Beta(θa,b)θa1(1θ)b1p(\theta) = \text{Beta}(\theta | a, b) \propto \theta^{a-1} (1-\theta)^{b-1}

Les paramètres aa et bb contrôlent la forme de l’a priori. Pour a=b=2a = b = 2, l’a priori favorise des valeurs de θ\theta proches de 0.5.

Le logarithme de l’a posteriori (vraisemblance plus a priori) est:

logp(θD)N1logθ+N0log(1θ)+(a1)logθ+(b1)log(1θ)\log p(\theta | \mathcal{D}) \propto N_1 \log \theta + N_0 \log(1-\theta) + (a-1) \log \theta + (b-1) \log(1-\theta)

En dérivant et en résolvant, l’estimateur MAP est:

θ^MAP=N1+a1N1+N0+a+b2\hat{\theta}_{\text{MAP}} = \frac{N_1 + a - 1}{N_1 + N_0 + a + b - 2}

Avec a=b=2a = b = 2 et nos 3 observations de faces:

θ^MAP=3+213+0+2+22=45=0.8\hat{\theta}_{\text{MAP}} = \frac{3 + 2 - 1}{3 + 0 + 2 + 2 - 2} = \frac{4}{5} = 0.8

Cette estimation est plus raisonnable que l’EMV θ^EMV=1\hat{\theta}_{\text{EMV}} = 1. L’a priori “tire” l’estimation vers des valeurs moins extrêmes.

Le choix a=b=2a = b = 2 correspond au lissage de Laplace (ou add-one smoothing): c’est comme si nous avions observé une face et une pile supplémentaires avant de commencer. Cette technique est particulièrement utile quand certains événements n’ont jamais été observés dans les données.

Source
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

theta_grid = np.linspace(0.001, 0.999, 200)

# Left: Prior, likelihood, posterior
ax = axes[0]
n1, n0 = 3, 0  # 3 heads, 0 tails
a, b = 2, 2    # Beta prior parameters

# Prior
prior = stats.beta.pdf(theta_grid, a, b)
prior = prior / prior.max()

# Likelihood
likelihood = theta_grid**n1 * (1 - theta_grid)**n0
likelihood = likelihood / likelihood.max()

# Posterior (Beta(a + n1, b + n0))
posterior = stats.beta.pdf(theta_grid, a + n1, b + n0)
posterior = posterior / posterior.max()

ax.plot(theta_grid, prior, 'g-', linewidth=2, label='A priori Beta(2,2)')
ax.plot(theta_grid, likelihood, 'b--', linewidth=2, label='Vraisemblance')
ax.plot(theta_grid, posterior, 'r-', linewidth=2, label='A posteriori')

theta_mle = n1 / (n1 + n0)
theta_map = (n1 + a - 1) / (n1 + n0 + a + b - 2)
ax.axvline(theta_mle, color='b', linestyle=':', alpha=0.7, label=f'EMV: {theta_mle:.2f}')
ax.axvline(theta_map, color='r', linestyle=':', alpha=0.7, label=f'MAP: {theta_map:.2f}')

ax.set_xlabel(r'$\theta$')
ax.set_ylabel('Densité (normalisée)')
ax.set_title('3 faces, 0 pile')
ax.legend(fontsize=8)
ax.set_xlim(0, 1)

# Right: Effect of different priors
ax = axes[1]
priors = [(1, 1, 'Uniforme'), (2, 2, 'Beta(2,2)'), (5, 5, 'Beta(5,5)')]

for a, b, label in priors:
    theta_map = (n1 + a - 1) / (n1 + n0 + a + b - 2)
    posterior = stats.beta.pdf(theta_grid, a + n1, b + n0)
    posterior = posterior / posterior.max()
    ax.plot(theta_grid, posterior, linewidth=2, label=f'{label}: MAP={theta_map:.2f}')

ax.axvline(1.0, color='gray', linestyle='--', alpha=0.5, label='EMV: 1.00')
ax.set_xlabel(r'$\theta$')
ax.set_ylabel('A posteriori (normalisé)')
ax.set_title('Effet de différents a priori')
ax.legend(fontsize=8)
ax.set_xlim(0, 1)

plt.tight_layout()
<Figure size 1000x400 with 2 Axes>

La figure de gauche montre comment l’a posteriori combine l’a priori et la vraisemblance. L’a priori Beta(2,2) “tire” l’estimation vers 0.5, résultant en un MAP de 0.8 au lieu de l’EMV de 1.0. La figure de droite montre l’effet de différents a priori: plus l’a priori est fort (variance faible), plus l’estimation est proche de 0.5.

Régression ridge = MAP avec a priori gaussien

Appliquons maintenant ce cadre bayésien à la régression linéaire. Si nous plaçons un a priori gaussien isotrope sur les paramètres:

p(θ)=N(θ0,σθ2I)p(\boldsymbol{\theta}) = \mathcal{N}(\boldsymbol{\theta} | \mathbf{0}, \sigma_\theta^2 \mathbf{I})

cet a priori exprime la croyance que les paramètres sont probablement proches de zéro, avec une incertitude contrôlée par σθ2\sigma_\theta^2.

Le logarithme négatif de cet a priori est:

logp(θ)=12σθ2θ22+constante-\log p(\boldsymbol{\theta}) = \frac{1}{2\sigma_\theta^2} \|\boldsymbol{\theta}\|_2^2 + \text{constante}

L’estimateur MAP devient:

θ^MAP=argminθ[NLV(θ)+12σθ2θ22]\hat{\boldsymbol{\theta}}_{\text{MAP}} = \arg\min_{\boldsymbol{\theta}} \left[ \text{NLV}(\boldsymbol{\theta}) + \frac{1}{2\sigma_\theta^2}\|\boldsymbol{\theta}\|_2^2 \right]

C’est exactement la régression ridge, avec λ=1/(2σθ2)\lambda = 1/(2\sigma_\theta^2). Cette correspondance nous donne une interprétation de l’hyperparamètre:

L’a priori gaussien sur les paramètres est parfois appelé dégradation des poids (weight decay) dans le contexte des réseaux de neurones, car il “tire” les paramètres vers zéro pendant l’entraînement.

Unification: deux langages pour un même problème

Les sections précédentes ont présenté deux approches pour l’apprentissage supervisé. La première, fondée sur la théorie de la décision, définit une fonction de perte et minimise le risque empirique. La seconde, probabiliste, modélise la distribution des données et estime les paramètres par maximum de vraisemblance ou maximum a posteriori.

Ces deux approches semblent différentes, mais elles aboutissent aux mêmes algorithmes. En choisissant la perte logarithmique (y,y^)=logp(yy^)\ell(y, \hat{y}) = -\log p(y | \hat{y}), le risque empirique devient exactement la log-vraisemblance négative (à un facteur 1/N1/N près). Minimiser l’un revient à minimiser l’autre. Sous bruit gaussien, cette perte se réduit à la perte quadratique; sous modèle de Bernoulli, à l’entropie croisée.

De même, ajouter une régularisation 2\ell_2 au risque empirique revient à supposer un a priori gaussien sur les paramètres. La régression ridge n’est rien d’autre que l’estimation MAP avec cet a priori. Le coefficient λ\lambda encode la force de notre croyance a priori: plus λ\lambda est grand, plus nous “tirons” les paramètres vers zéro.

Pourquoi alors utiliser deux langages? Parce qu’ils éclairent des aspects différents du problème. Le langage décisionnel (risque, perte, minimisation) est opérationnel: il dit comment construire un algorithme. Le langage probabiliste (vraisemblance, a priori, a posteriori) est interprétatif: il dit ce que nous supposons sur les données et pourquoi nos choix sont raisonnables. Ensemble, ils permettent de concevoir des algorithmes et de comprendre leur comportement.

Interprétation informationnelle

La théorie de l’information offre une troisième perspective. L’EMV peut se comprendre comme la recherche du modèle paramétrique le plus proche de la distribution empirique des données.

La distribution empirique place une masse 1/N1/N sur chaque observation:

pD(y)=1Ni=1Nδ(yyi)p_{\mathcal{D}}(y) = \frac{1}{N} \sum_{i=1}^N \delta(y - y_i)

La divergence de Kullback-Leibler mesure la dissimilarité entre deux distributions:

DKL(pq)=yp(y)logp(y)q(y)D_{\text{KL}}(p \| q) = \sum_y p(y) \log \frac{p(y)}{q(y)}

Cette quantité est toujours positive ou nulle, et vaut zéro si et seulement si les deux distributions sont identiques. En posant p=pDp = p_{\mathcal{D}} (ce que nous avons observé) et q=p(θ)q = p(\cdot | \boldsymbol{\theta}) (notre modèle), on peut montrer que:

argminθDKL(pDp(θ))=argminθNLV(θ)\arg\min_{\boldsymbol{\theta}} D_{\text{KL}}(p_{\mathcal{D}} \| p(\cdot|\boldsymbol{\theta})) = \arg\min_{\boldsymbol{\theta}} \text{NLV}(\boldsymbol{\theta})

L’EMV trouve les paramètres qui rendent notre modèle aussi proche que possible de ce que nous avons observé, au sens de la divergence KL. Cette interprétation géométrique complète les perspectives décisionnelle et probabiliste.

Résumé

Ce chapitre a développé le cadre formel de l’apprentissage supervisé à travers deux perspectives complémentaires.

La perspective décisionnelle définit le risque comme l’erreur moyenne attendue sur de nouvelles données, puis le remplace par le risque empirique (calculable sur les données d’entraînement). La minimisation du risque empirique cherche le modèle qui minimise cette erreur d’entraînement, en espérant qu’il généralisera. La régularisation pénalise la complexité pour éviter le surapprentissage.

La perspective probabiliste modélise explicitement comment les données ont été générées. Le maximum de vraisemblance (EMV) trouve les paramètres qui rendent les observations les plus probables; le maximum a posteriori (MAP) incorpore des croyances a priori sur les paramètres.

Ces deux perspectives aboutissent aux mêmes algorithmes. Sous bruit gaussien, l’EMV coïncide avec les moindres carrés. La régularisation 2\ell_2 correspond au MAP avec a priori gaussien sur les paramètres. Le langage décisionnel est opérationnel (comment construire l’algorithme); le langage probabiliste est interprétatif (pourquoi ces choix sont raisonnables).

Le chapitre suivant développe les outils théoriques pour quantifier quand et comment le risque empirique prédit le vrai risque.

Exercices

  1. Validation croisée:

    from sklearn.model_selection import KFold
    from sklearn.linear_model import Ridge
    
    lambdas = np.logspace(-4, 2, 20)
    kf = KFold(n_splits=5, shuffle=True, random_state=42)
    
    cv_scores = []
    for lam in lambdas:
        fold_scores = []
        for train_idx, val_idx in kf.split(X):
            model = Ridge(alpha=lam, fit_intercept=False)
            model.fit(X[train_idx], y[train_idx])
            y_pred = model.predict(X[val_idx])
            mse = np.mean((y[val_idx] - y_pred)**2)
            fold_scores.append(mse)
        cv_scores.append(np.mean(fold_scores))
  2. Visualisation et choix de λ:

    plt.plot(np.log10(lambdas), cv_scores)
    plt.xlabel('log10(λ)')
    plt.ylabel('EQM de validation')
    best_idx = np.argmin(cv_scores)
    best_lambda = lambdas[best_idx]

    La courbe montre typiquement:

    • EQM élevé pour λ\lambda très petit (surapprentissage)

    • EQM minimal pour λ\lambda intermédiaire

    • EQM qui remonte pour λ\lambda grand (sous-apprentissage)

    Le λ\lambda optimal se situe au minimum de la courbe (souvent autour de 10-1 à 100).

  3. Comparaison finale:

    MCO avec degré 10 surapprend fortement et a un EQM de test élevé. Ridge avec λ\lambda optimal a un EQM de test beaucoup plus faible car la régularisation empêche les coefficients d’exploser.


````{admonition} Exercice 14: Prédicteur de Bayes optimal (perte quadratique) ★★
:class: hint dropdown

Le prédicteur de Bayes optimal minimise le risque pour une perte donnée, en supposant que la vraie distribution $p(y|x)$ est connue.

Considérons la distribution conditionnelle suivante pour un $x$ donné:

$$
p(y|x) = 0.3 \cdot \mathcal{N}(y | 1, 0.5^2) + 0.7 \cdot \mathcal{N}(y | 4, 1^2)
$$

C'est un mélange de deux gaussiennes.

1. Tracez cette distribution $p(y|x)$.

2. Calculez l'espérance $\mathbb{E}[y|x]$ (utiliser la linéarité de l'espérance).

3. Pour la perte quadratique, le prédicteur optimal est la moyenne conditionnelle. Quelle est donc la prédiction optimale $\hat{y}^*$?

4. Calculez le risque de Bayes (l'erreur minimale atteignable): $\mathcal{R}^* = \mathbb{E}[(y - \hat{y}^*)^2 | x]$.

5. Si vous prédisiez le mode (la valeur la plus probable) au lieu de la moyenne, quelle serait votre prédiction? Quel serait le risque?
References
  1. Murphy, K. P. (2022). Probabilistic Machine Learning: An Introduction. MIT Press. https://probml.github.io/pml-book/book1.html